From ad160f121dd12328dcac0f69cdadc6ee3d45df28 Mon Sep 17 00:00:00 2001 From: Raffaele Montuoro Date: Fri, 29 Mar 2019 14:29:31 +0000 Subject: [PATCH] FV3: this commit #refs 61047 Merge FV3 chemistry development branch into master branch. Change-Id: I8014391ea97416ec8d757aa32ed38561b03aa679 --- atmos_model.F90 | 224 ++++++++++++++++++-- cpl/module_cap_cpl.F90 | 110 +++++++--- cpl/module_cplfields.F90 | 52 +++-- fv3_cap.F90 | 104 ++++++--- gfsphysics/GFS_layer/GFS_diagnostics.F90 | 157 +++++++++++++- gfsphysics/GFS_layer/GFS_physics_driver.F90 | 50 +++-- gfsphysics/GFS_layer/GFS_typedefs.F90 | 171 +++++++++++++-- gfsphysics/physics/gfdl_cloud_microphys.F90 | 8 +- module_fcst_grid_comp.F90 | 21 +- namphysics/NAM_layer/NAM_diagnostics.F90 | 156 +++++++++++++- namphysics/NAM_layer/NAM_typedefs.F90 | 153 +++++++++++-- namphysics/physics/gfdl_cloud_microphys.F90 | 8 +- 12 files changed, 1061 insertions(+), 153 deletions(-) diff --git a/atmos_model.F90 b/atmos_model.F90 index 423ccbb1b..59a8a6ea6 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -61,7 +61,7 @@ module atmos_model_mod operator(+), operator(-),real_to_time_type use field_manager_mod, only: MODEL_ATMOS use tracer_manager_mod, only: get_number_tracers, get_tracer_names, & - get_tracer_index + get_tracer_index, NO_TRACER use xgrid_mod, only: grid_box_type use atmosphere_mod, only: atmosphere_init use atmosphere_mod, only: atmosphere_restart @@ -792,20 +792,71 @@ end subroutine atmos_model_restart ! Retrieve ungridded dimensions of atmospheric model arrays ! -subroutine get_atmos_model_ungridded_dim(nlev, ntracers, nsoillev) +subroutine get_atmos_model_ungridded_dim(nlev, nsoillev, ntracers, & + num_diag_sfc_emis_flux, num_diag_down_flux, num_diag_type_down_flux, & + num_diag_burn_emis_flux, num_diag_cmass) - integer, optional, intent(out) :: nlev, ntracers, nsoillev + integer, optional, intent(out) :: nlev, nsoillev, ntracers, & + num_diag_sfc_emis_flux, num_diag_down_flux, num_diag_type_down_flux, & + num_diag_burn_emis_flux, num_diag_cmass - if (present(nlev)) nlev = Atm_block%npz - if (present(ntracers)) call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) + !--- number of atmospheric vertical levels + if (present(nlev)) nlev = Atm_block%npz + + !--- number of soil levels if (present(nsoillev)) then nsoillev = 0 if (allocated(IPD_Data)) then if (associated(IPD_Data(1)%Sfcprop%slc)) & - nsoillev = size(IPD_Data(1)%Sfcprop%slc, 2) + nsoillev = size(IPD_Data(1)%Sfcprop%slc, dim=2) end if end if + !--- total number of atmospheric tracers + if (present(ntracers)) call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) + + !--- number of tracers used in chemistry diagnostic output + if (present(num_diag_down_flux)) then + num_diag_down_flux = 0 + if (associated(IPD_Data(1)%IntDiag%sedim)) & + num_diag_down_flux = size(IPD_Data(1)%IntDiag%sedim, dim=2) + if (present(num_diag_type_down_flux)) then + num_diag_type_down_flux = 0 + if (associated(IPD_Data(1)%IntDiag%sedim)) & + num_diag_type_down_flux = num_diag_type_down_flux + 1 + if (associated(IPD_Data(1)%IntDiag%drydep)) & + num_diag_type_down_flux = num_diag_type_down_flux + 1 + if (associated(IPD_Data(1)%IntDiag%wetdpl)) & + num_diag_type_down_flux = num_diag_type_down_flux + 1 + if (associated(IPD_Data(1)%IntDiag%wetdpc)) & + num_diag_type_down_flux = num_diag_type_down_flux + 1 + end if + end if + + !--- number of bins for chemistry diagnostic output + if (present(num_diag_sfc_emis_flux)) then + num_diag_sfc_emis_flux = 0 + if (associated(IPD_Data(1)%IntDiag%duem)) & + num_diag_sfc_emis_flux = size(IPD_Data(1)%IntDiag%duem, dim=2) + if (associated(IPD_Data(1)%IntDiag%ssem)) & + num_diag_sfc_emis_flux = & + num_diag_sfc_emis_flux + size(IPD_Data(1)%IntDiag%ssem, dim=2) + end if + + !--- number of tracers used in emission diagnostic output + if (present(num_diag_burn_emis_flux)) then + num_diag_burn_emis_flux = 0 + if (associated(IPD_Data(1)%IntDiag%abem)) & + num_diag_burn_emis_flux = size(IPD_Data(1)%IntDiag%abem, dim=2) + end if + + !--- number of tracers used in column mass density diagnostics + if (present(num_diag_cmass)) then + num_diag_cmass = 0 + if (associated(IPD_Data(1)%IntDiag%aecm)) & + num_diag_cmass = size(IPD_Data(1)%IntDiag%aecm, dim=2) + end if + end subroutine get_atmos_model_ungridded_dim ! @@ -836,21 +887,22 @@ subroutine update_atmos_chemistry(state, rc) !--- local variables integer :: localrc - integer :: ni, nj, nk, nt, ntoz + integer :: ni, nj, nk, nt, ntb, nte integer :: nb, ix, i, j, k, it integer :: ib, jb - real(ESMF_KIND_R8), dimension(:,:,:), pointer :: prsl, phil, & - prsi, phii, & - temp, & + real(ESMF_KIND_R8), dimension(:,:,:), pointer :: prsl, phil, & + prsi, phii, & + temp, & ua, va, vvl, & - dkt, slc - real(ESMF_KIND_R8), dimension(:,:,:,:), pointer :: q + dkt, slc, & + qb, qm, qu + real(ESMF_KIND_R8), dimension(:,:,:,:), pointer :: qd, q real(ESMF_KIND_R8), dimension(:,:), pointer :: hpbl, area, stype, rainc, & uustar, rain, sfcdsw, slmsk, tsfc, shfsfc, snowd, vtype, vfrac, zorl - logical, parameter :: diag = .true. +! logical, parameter :: diag = .true. ! -- begin if (present(rc)) rc = ESMF_SUCCESS @@ -858,6 +910,8 @@ subroutine update_atmos_chemistry(state, rc) ni = Atm_block%iec - Atm_block%isc + 1 nj = Atm_block%jec - Atm_block%jsc + 1 nk = Atm_block%npz + + !--- get total number of tracers call get_number_tracers(MODEL_ATMOS, num_tracers=nt) select case (trim(state)) @@ -867,12 +921,37 @@ subroutine update_atmos_chemistry(state, rc) farrayPtr4d=q, rc=localrc) if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__, rcToReturn=rc)) return + call cplFieldGet(state,'inst_tracer_up_surface_flx', & + farrayPtr3d=qu, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, rcToReturn=rc)) return + call cplFieldGet(state,'inst_tracer_down_surface_flx', & + farrayPtr4d=qd, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, rcToReturn=rc)) return + call cplFieldGet(state,'inst_tracer_clmn_mass_dens', & + farrayPtr3d=qm, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, rcToReturn=rc)) return + call cplFieldGet(state,'inst_tracer_anth_biom_flx', & + farrayPtr3d=qb, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, rcToReturn=rc)) return - !--- tracers quantities - !--- locate the end location of standard atmospheric tracers, marked by ozone - ntoz = get_tracer_index(MODEL_ATMOS, 'o3mr') + !--- do not import tracer concentrations by default + ntb = nt + 1 + nte = nt + + !--- if chemical tracers are present, set bounds appropriately + if (IPD_Control%ntchm > 0) then + if (IPD_Control%ntchs /= NO_TRACER) then + ntb = IPD_Control%ntchs + nte = IPD_Control%ntchm + ntb - 1 + end if + end if - do it = ntoz + 1, nt + !--- tracer concentrations + do it = ntb, nte !$OMP parallel do default (none) & !$OMP shared (it, nk, nj, ni, Atm_block, IPD_Data, q) & !$OMP private (k, j, jb, i, ib, nb, ix) @@ -889,9 +968,97 @@ subroutine update_atmos_chemistry(state, rc) enddo enddo - if (diag) then + !--- tracer diagnostics + !--- (a) column mass densities + do it = 1, size(qm, dim=3) +!$OMP parallel do default (none) & +!$OMP shared (it, nj, ni, Atm_block, IPD_Data, qm) & +!$OMP private (j, jb, i, ib, nb, ix) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + IPD_Data(nb)%IntDiag%aecm(ix,it) = qm(i,j,it) + enddo + enddo + enddo + + !--- (b) dust and sea salt emissions + ntb = size(IPD_Data(1)%IntDiag%duem, dim=2) + nte = size(qu, dim=3) + do it = 1, min(ntb, nte) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + IPD_Data(nb)%IntDiag%duem(ix,it) = qu(i,j,it) + enddo + enddo + enddo + + nte = nte - ntb + do it = 1, min(size(IPD_Data(1)%IntDiag%ssem, dim=2), nte) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + IPD_Data(nb)%IntDiag%ssem(ix,it) = qu(i,j,it+ntb) + enddo + enddo + enddo + + !--- (c) sedimentation and dry/wet deposition + do it = 1, size(qd, dim=3) +!$OMP parallel do default (none) & +!$OMP shared (it, nj, ni, Atm_block, IPD_Data, qd) & +!$OMP private (j, jb, i, ib, nb, ix) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + IPD_Data(nb)%IntDiag%sedim (ix,it) = qd(i,j,it,1) + IPD_Data(nb)%IntDiag%drydep(ix,it) = qd(i,j,it,2) + IPD_Data(nb)%IntDiag%wetdpl(ix,it) = qd(i,j,it,3) + IPD_Data(nb)%IntDiag%wetdpc(ix,it) = qd(i,j,it,4) + enddo + enddo + enddo + + !--- (d) anthropogenic and biomass burning emissions + do it = 1, size(qb, dim=3) +!$OMP parallel do default (none) & +!$OMP shared (it, nj, ni, Atm_block, IPD_Data, qb) & +!$OMP private (j, jb, i, ib, nb, ix) + do j = 1, nj + jb = j + Atm_block%jsc - 1 + do i = 1, ni + ib = i + Atm_block%isc - 1 + nb = Atm_block%blkno(ib,jb) + ix = Atm_block%ixp(ib,jb) + IPD_Data(nb)%IntDiag%abem(ix,it) = qb(i,j,it) + enddo + enddo + enddo + + if (IPD_Control%debug) then write(6,'("update_atmos: ",a,": qgrs - min/max/avg",3g16.6)') & trim(state), minval(q), maxval(q), sum(q)/size(q) + write(6,'("update_atmos: ",a,": qup - min/max/avg",3g16.6)') & + trim(state), minval(qu), maxval(qu), sum(qu)/size(qu) + write(6,'("update_atmos: ",a,": qdwn - min/max/avg",3g16.6)') & + trim(state), minval(qd), maxval(qd), sum(qd)/size(qd) + write(6,'("update_atmos: ",a,": qcmd - min/max/avg",3g16.6)') & + trim(state), minval(qm), maxval(qm), sum(qm)/size(qm) + write(6,'("update_atmos: ",a,": qabb - min/max/avg",3g16.6)') & + trim(state), minval(qb), maxval(qb), sum(qb)/size(qb) end if case ('export') @@ -935,7 +1102,8 @@ subroutine update_atmos_chemistry(state, rc) if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__, rcToReturn=rc)) return - call cplFieldGet(state,'inst_soil_moisture_content', farrayPtr3d=slc, rc=rc) + call cplFieldGet(state,'inst_soil_moisture_content', & + farrayPtr3d=slc, rc=localrc) if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__, rcToReturn=rc)) return @@ -1078,7 +1246,8 @@ subroutine update_atmos_chemistry(state, rc) area(i,j) = IPD_Data(nb)%Grid%area(ix) stype(i,j) = IPD_Data(nb)%Sfcprop%stype(ix) rainc(i,j) = IPD_Data(nb)%Coupling%rainc_cpl(ix) - rain(i,j) = IPD_Data(nb)%Coupling%rain_cpl(ix) + rain(i,j) = IPD_Data(nb)%Coupling%rain_cpl(ix) & + + IPD_Data(nb)%Coupling%snow_cpl(ix) uustar(i,j) = IPD_Data(nb)%Sfcprop%uustar(ix) sfcdsw(i,j) = IPD_Data(nb)%Coupling%sfcdsw(ix) slmsk(i,j) = IPD_Data(nb)%Sfcprop%slmsk(ix) @@ -1092,7 +1261,7 @@ subroutine update_atmos_chemistry(state, rc) enddo enddo - if (diag) then + if (IPD_Control%debug) then ! -- diagnostics write(6,'("update_atmos: prsi - min/max/avg",3g16.6)') minval(prsi), maxval(prsi), sum(prsi)/size(prsi) write(6,'("update_atmos: phii - min/max/avg",3g16.6)') minval(phii), maxval(phii), sum(phii)/size(phii) @@ -1208,7 +1377,7 @@ subroutine assign_importdata(rc) real(kind=ESMF_KIND_R4), dimension(:,:), pointer :: datar42d real(kind=ESMF_KIND_R8), dimension(:,:), pointer :: datar82d real(kind=IPD_kind_phys), dimension(:,:), pointer :: datar8 - logical found, lcpl_fice + logical found, isFieldCreated, lcpl_fice ! !------------------------------------------------------------------------------ ! @@ -1230,16 +1399,25 @@ subroutine assign_importdata(rc) ! Each import field is only available if it was connected in the ! import state. found = .false. - if (ESMF_FieldIsCreated(importFields(n))) then + + isFieldCreated = ESMF_FieldIsCreated(importFields(n), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + + if (isFieldCreated) then ! put the data from local cubed sphere grid to column grid for phys datar8 = -99999.0 call ESMF_FieldGet(importFields(n), dimCount=dimCount ,typekind=datatype, & name=impfield_name, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out if ( dimCount == 2) then if ( datatype == ESMF_TYPEKIND_R8) then call ESMF_FieldGet(importFields(n),farrayPtr=datar82d,localDE=0, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out datar8=datar82d if (mpp_pe() == mpp_root_pe() .and. debug) print *,'in cplIMP,atmos gets ',trim(impfield_name),' datar8=',maxval(datar8),minval(datar8), & datar8(isc,jsc) diff --git a/cpl/module_cap_cpl.F90 b/cpl/module_cap_cpl.F90 index 9855fc2d1..3f9257765 100644 --- a/cpl/module_cap_cpl.F90 +++ b/cpl/module_cap_cpl.F90 @@ -71,7 +71,8 @@ subroutine realizeConnectedInternCplField(state, field, standardName, grid, rc) ! local variables character(len=80) :: fieldName type(ESMF_ArraySpec) :: arrayspec - integer :: i + integer :: i, localrc + logical :: isConnected real(ESMF_KIND_R8), pointer :: fptr(:,:) if (present(rc)) rc = ESMF_SUCCESS @@ -79,53 +80,66 @@ subroutine realizeConnectedInternCplField(state, field, standardName, grid, rc) fieldName = standardName ! use standard name as field name !! Create fields using wam2dmesh if they are WAM fields - if (NUOPC_IsConnected(state, fieldName=fieldName)) then + isConnected = NUOPC_IsConnected(state, fieldName=fieldName, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out - field = ESMF_FieldCreate(grid, ESMF_TYPEKIND_R8, name=fieldName, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + if (isConnected) then + + field = ESMF_FieldCreate(grid, ESMF_TYPEKIND_R8, name=fieldName, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & - file=__FILE__)) & - return ! bail out - call NUOPC_Realize(state, field=field, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + call NUOPC_Realize(state, field=field, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & - file=__FILE__)) & - return + file=__FILE__, & + rcToReturn=rc)) return - call ESMF_FieldGet(field, farrayPtr=fptr, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + call ESMF_FieldGet(field, farrayPtr=fptr, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & - file=__FILE__)) & - return ! bail out + file=__FILE__, & + rcToReturn=rc)) return ! bail out fptr=0.d0 ! zero out the entire field - call NUOPC_SetAttribute(field, name="Updated", value="true", rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + call NUOPC_SetAttribute(field, name="Updated", value="true", rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & - file=__FILE__)) & - return ! bail out + file=__FILE__, & + rcToReturn=rc)) return ! bail out else ! remove a not connected Field from State - call ESMF_StateRemove(state, (/fieldName/), rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + call ESMF_StateRemove(state, (/fieldName/), rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & - file=__FILE__)) & - return ! bail out + file=__FILE__, & + rcToReturn=rc)) return ! bail out endif end subroutine realizeConnectedInternCplField !----------------------------------------------------------------------------- - subroutine realizeConnectedCplFields(state, grid, & - numLevels, numSoilLayers, numTracers, fieldNames, fieldTypes, fieldList, rc) + subroutine realizeConnectedCplFields(state, grid, & + numLevels, numSoilLayers, numTracers, num_diag_sfc_emis_flux, & + num_diag_down_flux, num_diag_type_down_flux, num_diag_burn_emis_flux, & + num_diag_cmass, fieldNames, fieldTypes, fieldList, rc) type(ESMF_State), intent(inout) :: state type(ESMF_Grid), intent(in) :: grid integer, intent(in) :: numLevels integer, intent(in) :: numSoilLayers integer, intent(in) :: numTracers + integer, intent(in) :: num_diag_sfc_emis_flux + integer, intent(in) :: num_diag_down_flux + integer, intent(in) :: num_diag_type_down_flux + integer, intent(in) :: num_diag_burn_emis_flux + integer, intent(in) :: num_diag_cmass character(len=*), dimension(:), intent(in) :: fieldNames character(len=*), dimension(:), intent(in) :: fieldTypes type(ESMF_Field), dimension(:), intent(out) :: fieldList @@ -133,8 +147,8 @@ subroutine realizeConnectedCplFields(state, grid, & ! local variables integer :: item + logical :: isConnected type(ESMF_Field) :: field - real(ESMF_KIND_R8), pointer :: fptr(:,:) ! begin rc = ESMF_SUCCESS @@ -147,7 +161,10 @@ subroutine realizeConnectedCplFields(state, grid, & end if do item = 1, size(fieldNames) - if (NUOPC_IsConnected(state, fieldName=trim(fieldNames(item)))) then + isConnected = NUOPC_IsConnected(state, fieldName=trim(fieldNames(item)), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + if (isConnected) then select case (fieldTypes(item)) case ('l','layer') field = ESMF_FieldCreate(grid, typekind=ESMF_TYPEKIND_R8, & @@ -164,7 +181,32 @@ subroutine realizeConnectedCplFields(state, grid, & case ('t','tracer') field = ESMF_FieldCreate(grid, typekind=ESMF_TYPEKIND_R8, & name=trim(fieldNames(item)), & - ungriddedLBound=(/1,1/), ungriddedUBound=(/numLevels, numTracers/), rc=rc) + ungriddedLBound=(/1, 1/), ungriddedUBound=(/numLevels, numTracers/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + case ('u','tracer_up_flux') + field = ESMF_FieldCreate(grid, typekind=ESMF_TYPEKIND_R8, & + name=trim(fieldNames(item)), & + ungriddedLBound=(/1/), ungriddedUBound=(/num_diag_sfc_emis_flux/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + case ('d','tracer_down_flx') + field = ESMF_FieldCreate(grid, typekind=ESMF_TYPEKIND_R8, & + name=trim(fieldNames(item)), & + ungriddedLBound=(/1, 1/), & + ungriddedUBound=(/num_diag_down_flux, num_diag_type_down_flux/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + case ('b','tracer_anth_biom_emission') + field = ESMF_FieldCreate(grid, typekind=ESMF_TYPEKIND_R8, & + name=trim(fieldNames(item)), & + ungriddedLBound=(/1/), ungriddedUBound=(/num_diag_burn_emis_flux/), rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + case ('c','tracer_column_mass_density') + field = ESMF_FieldCreate(grid, typekind=ESMF_TYPEKIND_R8, & + name=trim(fieldNames(item)), & + ungriddedLBound=(/1/), ungriddedUBound=(/num_diag_cmass/), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return ! bail out case ('s','surface') @@ -190,12 +232,11 @@ subroutine realizeConnectedCplFields(state, grid, & file=__FILE__)) & return ! -- zero out field - call ESMF_FieldGet(field, farrayPtr=fptr, rc=rc) + call ESMF_FieldFill(field, dataFillScheme="const", const1=0._ESMF_KIND_R8, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & return ! bail out - fptr=0._ESMF_KIND_R8 ! zero out the entire field ! -- save field fieldList(item) = field else @@ -344,7 +385,8 @@ subroutine ESMFPP_RegridWrite(inField, outGrid, regridMethod, fileName, fieldNam integer, intent(in) :: timeslice integer, intent(inout) :: rc - ! local arguments + ! local variables + integer :: srcTermProcessing type(ESMF_Routehandle) :: rh type(ESMF_Field) :: outField @@ -354,10 +396,13 @@ subroutine ESMFPP_RegridWrite(inField, outGrid, regridMethod, fileName, fieldNam file=__FILE__)) & return ! bail out + ! Perform entire regridding arithmetic on the destination PET + srcTermProcessing = 0 ! For other options for the regrid operation, please refer to: ! http://www.earthsystemmodeling.org/esmf_releases/last_built/ESMF_refdoc/node5.html#SECTION050366000000000000000 call ESMF_FieldRegridStore(inField, outField, regridMethod=regridMethod, & unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, & + srcTermProcessing=srcTermProcessing, & Routehandle=rh, & rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & @@ -365,7 +410,10 @@ subroutine ESMFPP_RegridWrite(inField, outGrid, regridMethod, fileName, fieldNam file=__FILE__)) & return ! bail out - call ESMF_FieldRegrid(inField, outField, Routehandle=rh, rc=rc) + ! Use fixed ascending order for the sum terms based on their source + ! sequence index to ensure bit-for-bit reproducibility + call ESMF_FieldRegrid(inField, outField, Routehandle=rh, & + termorderflag=ESMF_TERMORDER_SRCSEQ, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & diff --git a/cpl/module_cplfields.F90 b/cpl/module_cplfields.F90 index 216177776..38de36871 100644 --- a/cpl/module_cplfields.F90 +++ b/cpl/module_cplfields.F90 @@ -117,19 +117,19 @@ module module_cplfields ! Set exportFieldShare to .true. if field is provided as memory reference ! to coupled components logical, public, parameter :: exportFieldShare(NexportFields) = (/ & - .true., .true., .true., .true., .true., & - .true., .true., .true., .true., .true., & - .true., .true., .true., .true., .true., & - .true., .true., .true., .true., .true., & - .true., .true., .true., .true., .true., & - .false.,.false.,.false.,.false.,.false., & - .false.,.false.,.false.,.false.,.false., & - .false.,.false.,.false.,.false.,.false., & + .true. ,.true. ,.true. ,.true. ,.true. , & + .true. ,.true. ,.true. ,.true. ,.true. , & + .true. ,.true. ,.true. ,.true. ,.true. , & + .true. ,.true. ,.true. ,.true. ,.true. , & + .true. ,.true. ,.false.,.false.,.false., & .false.,.false.,.false.,.false.,.false., & + .false.,.false.,.false.,.false.,.true. , & + .false.,.false.,.false.,.false.,.true. , & .false.,.false.,.false.,.false.,.false., & .false.,.false.,.false.,.false.,.false., & .false.,.false.,.false.,.false.,.false., & .false.,.false.,.false.,.false.,.false., & + .false.,.false.,.true. ,.false.,.false., & .false.,.false.,.false.,.false.,.false. & ! .false.,.false.,.false.,.false.,.false., & ! .false.,.false.,.false. & @@ -137,11 +137,15 @@ module module_cplfields real(kind=8), allocatable, public :: exportData(:,:,:) ! Import Fields ---------------------------------------- - integer, public, parameter :: NimportFields = 12 + integer, public, parameter :: NimportFields = 16 logical, public :: importFieldsValid(NimportFields) type(ESMF_Field), target, public :: importFields(NimportFields) character(len=*), public, parameter :: importFieldsList(NimportFields) = (/ & "inst_tracer_mass_frac ", & + "inst_tracer_up_surface_flx ", & + "inst_tracer_down_surface_flx ", & + "inst_tracer_clmn_mass_dens ", & + "inst_tracer_anth_biom_flx ", & "land_mask ", & "sea_ice_surface_temperature ", & "sea_surface_temperature ", & @@ -160,18 +164,18 @@ module module_cplfields "mean_snow_volume " & /) character(len=*), public, parameter :: importFieldTypes(NimportFields) = (/ & - "t", & - "s","s","s","s", & - "s","s","s", & - "s","s","s","s" & + "t","u","d","c","b", & + "s","s","s","s","s", & + "s","s","s","s","s", & + "s" & /) ! Set importFieldShare to .true. if field is provided as memory reference ! from coupled components logical, public, parameter :: importFieldShare(NimportFields) = (/ & - .true., & - .false.,.false.,.false.,.false., & - .false.,.false.,.false., & - .false.,.false.,.false.,.false. & + .true. ,.true. ,.true. ,.true. ,.true. , & + .false.,.false.,.false.,.false.,.false., & + .false.,.false.,.false.,.false.,.false., & + .false. & /) ! Methods @@ -190,6 +194,7 @@ subroutine fillExportFields(data_a2oi, rc) integer :: localrc integer :: n,dimCount + logical :: isCreated type(ESMF_TypeKind_Flag) :: datatype real(kind=ESMF_KIND_R4), dimension(:,:), pointer :: datar42d real(kind=ESMF_KIND_R8), dimension(:,:), pointer :: datar82d @@ -198,7 +203,11 @@ subroutine fillExportFields(data_a2oi, rc) if (present(rc)) rc=ESMF_SUCCESS do n=1, size(exportFields) - if (ESMF_FieldIsCreated(exportFields(n))) then + isCreated = ESMF_FieldIsCreated(exportFields(n), rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, & + rcToReturn=rc)) return + if (isCreated) then ! set data call ESMF_FieldGet(exportFields(n), dimCount=dimCount, & typekind=datatype, rc=localrc) @@ -312,6 +321,7 @@ subroutine cplFieldGet(state, name, localDe, & !--- local variables integer :: localrc integer :: de, item, fieldCount, rank + logical :: isCreated type(ESMF_Field), pointer :: fieldList(:) character(len=ESMF_MAXSTR) :: fieldName @@ -331,7 +341,11 @@ subroutine cplFieldGet(state, name, localDe, & rcToReturn=rc)) return do item = 1, fieldCount - if (NUOPC_IsConnected(fieldList(item))) then + isCreated = ESMF_FieldIsCreated(fieldList(item), rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__, & + rcToReturn=rc)) return + if (isCreated) then call ESMF_FieldGet(fieldList(item), name=fieldName, rc=localrc) if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__, & diff --git a/fv3_cap.F90 b/fv3_cap.F90 index 6258a4d25..80da17255 100644 --- a/fv3_cap.F90 +++ b/fv3_cap.F90 @@ -45,8 +45,12 @@ module fv3gfs_cap_mod stdlat1, stdlat2, dx, dy ! use module_fcst_grid_comp, only: fcstSS => SetServices, & - fcstGrid, numLevels, numTracers, & - numSoilLayers + fcstGrid, numLevels, numSoilLayers, & + numTracers, num_diag_sfc_emis_flux, & + num_diag_down_flux, & + num_diag_type_down_flux, & + num_diag_burn_emis_flux, num_diag_cmass + use module_wrt_grid_comp, only: wrtSS => SetServices ! use module_cplfields, only: nExportFields, exportFields, & @@ -138,8 +142,10 @@ subroutine SetServices(gcomp, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return ! bail out + ! temporarily disable timestamp check for import fields as chemistry fields + ! require to be handled differently due to the coupling run sequence call NUOPC_CompSpecialize(gcomp, specLabel=model_label_CheckImport, & - specRoutine=fv3_checkimport, rc=rc) + specRoutine=NUOPC_NoOp, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return ! bail out @@ -243,6 +249,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) integer :: mpi_comm_atm integer :: i, j, k, io_unit, urc integer :: petcount, mype + logical :: isPetLocal logical :: OPENED character(ESMF_MAXSTR) :: name logical :: fcstpe @@ -825,7 +832,12 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! --- advertise Fields in importState and exportState ------------------- if( cpl ) then - if (ESMF_GridCompIsPetLocal(fcstComp, rc=rc)) then + + isPetLocal = ESMF_GridCompIsPetLocal(fcstComp, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + + if (isPetLocal) then ! importable fields: do i = 1, size(ImportFieldsList) @@ -875,6 +887,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) integer, intent(out) :: rc ! local variables + logical :: isPetLocal integer :: n rc = ESMF_SUCCESS @@ -882,19 +895,26 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) ! --- conditionally realize or remove Fields in importState and exportState ------------------- if ( cpl ) then - if (ESMF_GridCompIsPetLocal(fcstComp, rc=rc)) then + + isPetLocal = ESMF_GridCompIsPetLocal(fcstComp, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + + if (isPetLocal) then ! -- realize connected fields in exportState - call realizeConnectedCplFields(exportState, fcstGrid, & - numLevels, numSoilLayers, numTracers, & - exportFieldsList, exportFieldTypes, exportFields, rc) + call realizeConnectedCplFields(exportState, fcstGrid, & + numLevels, numSoilLayers, numTracers, num_diag_sfc_emis_flux, & + num_diag_down_flux, num_diag_type_down_flux, num_diag_burn_emis_flux, & + num_diag_cmass, exportFieldsList, exportFieldTypes, exportFields, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return ! bail out ! -- realize connected fields in importState - call realizeConnectedCplFields(importState, fcstGrid, & - numLevels, numSoilLayers, numTracers, & - importFieldsList, importFieldTypes, importFields, rc) + call realizeConnectedCplFields(importState, fcstGrid, & + numLevels, numSoilLayers, numTracers, num_diag_sfc_emis_flux, & + num_diag_down_flux, num_diag_type_down_flux, num_diag_burn_emis_flux, & + num_diag_cmass, importFieldsList, importFieldTypes, importFields, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return ! bail out end if @@ -918,7 +938,7 @@ subroutine ModelAdvance(gcomp, rc) integer(ESMF_KIND_I8) :: n_interval, time_elapsed_sec ! integer :: na, i, urc - logical :: lalarm, reconcileFlag + logical :: isAlarmEnabled, isAlarmRinging, lalarm, reconcileFlag character(len=*),parameter :: subname='(fv3_cap:ModelAdvance)' character(240) :: msgString !jw debug @@ -1062,18 +1082,36 @@ subroutine ModelAdvance(gcomp, rc) if (nfhmax_hf > 0) then if(currtime <= starttime+output_hfmax) then - if(ESMF_AlarmIsEnabled(alarm = ALARM_OUTPUT_HF, rc = RC)) then - if( ESMF_AlarmIsRinging(alarm = ALARM_OUTPUT_HF,rc = Rc)) LALARM = .true. + isAlarmEnabled = ESMF_AlarmIsEnabled(alarm = ALARM_OUTPUT_HF, rc = RC) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + if(isAlarmEnabled) then + isAlarmRinging = ESMF_AlarmIsRinging(alarm = ALARM_OUTPUT_HF,rc = Rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + if (isAlarmRinging) LALARM = .true. endif else - if(ESMF_AlarmIsEnabled(alarm = ALARM_OUTPUT, rc = RC)) then - if(ESMF_AlarmIsRinging(alarm = ALARM_OUTPUT,rc = Rc)) LALARM = .true. + isAlarmEnabled = ESMF_AlarmIsEnabled(alarm = ALARM_OUTPUT, rc = RC) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + if(isAlarmEnabled) then + isAlarmRinging = ESMF_AlarmIsRinging(alarm = ALARM_OUTPUT,rc = Rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + if (isAlarmRinging) LALARM = .true. endif endif endif ! - if(ESMF_AlarmIsEnabled(alarm = ALARM_OUTPUT, rc = RC)) then - if(ESMF_AlarmIsRinging(alarm = ALARM_OUTPUT,rc = Rc)) LALARM = .true. + isAlarmEnabled = ESMF_AlarmIsEnabled(alarm = ALARM_OUTPUT, rc = RC) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + if(isAlarmEnabled) then + isAlarmRinging = ESMF_AlarmIsRinging(alarm = ALARM_OUTPUT,rc = Rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + if (isAlarmRinging) LALARM = .true. endif ! if (mype == 0 .or. mype == lead_wrttask(1)) print *,' aft fcst run lalarm=',lalarm, & ! 'FBcount=',FBcount,'na=',na @@ -1295,7 +1333,7 @@ subroutine ModelAdvance_phase2(gcomp, rc) integer(ESMF_KIND_I8) :: n_interval, time_elapsed_sec ! integer :: na, i, urc - logical :: lalarm, reconcileFlag + logical :: isAlarmEnabled, isAlarmRinging, lalarm, reconcileFlag character(len=*),parameter :: subname='(fv3_cap:ModelAdvance_phase2)' character(240) :: msgString !jw debug @@ -1364,19 +1402,35 @@ subroutine ModelAdvance_phase2(gcomp, rc) if (nfhmax_hf > 0) then if(currtime <= starttime+output_hfmax) then - if(ESMF_AlarmIsEnabled(alarm = ALARM_OUTPUT_HF, rc = RC)) then - if( ESMF_AlarmIsRinging(alarm = ALARM_OUTPUT_HF,rc = Rc)) LALARM = .true. + isAlarmEnabled = ESMF_AlarmIsEnabled(alarm = ALARM_OUTPUT_HF, rc = RC) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + if(isAlarmEnabled) then + isAlarmRinging = ESMF_AlarmIsRinging(alarm = ALARM_OUTPUT_HF,rc = Rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + if (isAlarmRinging) LALARM = .true. endif else - if(ESMF_AlarmIsEnabled(alarm = ALARM_OUTPUT, rc = RC)) then - if(ESMF_AlarmIsRinging(alarm = ALARM_OUTPUT,rc = Rc)) LALARM = .true. + isAlarmEnabled = ESMF_AlarmIsEnabled(alarm = ALARM_OUTPUT, rc = RC) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + if(isAlarmEnabled) then + isAlarmRinging = ESMF_AlarmIsRinging(alarm = ALARM_OUTPUT,rc = Rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + if (isAlarmRinging) LALARM = .true. endif endif endif ! - if(ESMF_AlarmIsEnabled(alarm = ALARM_OUTPUT, rc = RC)) then - if(ESMF_AlarmIsRinging(alarm = ALARM_OUTPUT,rc = Rc)) LALARM = .true. + isAlarmEnabled = ESMF_AlarmIsEnabled(alarm = ALARM_OUTPUT, rc = RC) + if(isAlarmEnabled) then + isAlarmRinging = ESMF_AlarmIsRinging(alarm = ALARM_OUTPUT,rc = Rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, file=__FILE__)) return ! bail out + if (isAlarmRinging) LALARM = .true. endif if (mype == 0 .or. mype == lead_wrttask(1)) print *,' aft fcst run lalarm=',lalarm, & 'FBcount=',FBcount,'na=',na diff --git a/gfsphysics/GFS_layer/GFS_diagnostics.F90 b/gfsphysics/GFS_layer/GFS_diagnostics.F90 index 0e57d281f..39e2491e1 100644 --- a/gfsphysics/GFS_layer/GFS_diagnostics.F90 +++ b/gfsphysics/GFS_layer/GFS_diagnostics.F90 @@ -98,7 +98,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop type(GFS_init_type), intent(in) :: Init_parm !--- local variables - integer :: idx, num, nb, nblks, NFXR + integer :: idt, idx, num, nb, nblks, NFXR character(len=2) :: xtra real(kind=kind_phys), parameter :: cn_one = 1._kind_phys real(kind=kind_phys), parameter :: cn_100 = 100._kind_phys @@ -2864,6 +2864,161 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop endif ! print *,'in gfdl_diag_register,af all extdiag, idx=',idx +! -- chemistry diagnostic variables + if (Model%cplchm) then + + if (Model%ntchm > 0) then + + if (associated(IntDiag(1)%duem)) then + do num = 1, size(IntDiag(1)%duem, dim=2) + idx = idx + 1 + ExtDiag(idx)%axes = 2 + write(ExtDiag(idx)%name,'("duem",i3.3)') num + write(ExtDiag(idx)%desc,'("Dust Emission Bin ",i0)') num + ExtDiag(idx)%unit = 'kg/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%duem(:,num) + enddo + enddo + endif + + if (associated(IntDiag(1)%ssem)) then + do num = 1, size(IntDiag(1)%ssem, dim=2) + idx = idx + 1 + ExtDiag(idx)%axes = 2 + write(ExtDiag(idx)%name,'("ssem",i3.3)') num + write(ExtDiag(idx)%desc,'("Seasalt Emission Bin ",i0)') num + ExtDiag(idx)%unit = 'kg/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%ssem(:,num) + enddo + enddo + endif + + if (associated(Model%ntdiag)) then + idt = 0 + do num = Model%ntchs, Model%ntchm + Model%ntchs - 1 + if (Model%ntdiag(num-Model%ntchs+1)) then + idt = idt + 1 + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = trim(Model%tracer_names(num)) // 'sd' + ExtDiag(idx)%desc = trim(Model%tracer_names(num)) // ' Sedimentation' + ExtDiag(idx)%unit = 'kg/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%sedim(:,idt) + enddo + + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = trim(Model%tracer_names(num)) // 'dp' + ExtDiag(idx)%desc = trim(Model%tracer_names(num)) // ' Dry Deposition' + ExtDiag(idx)%unit = 'kg/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%drydep(:,idt) + enddo + + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = trim(Model%tracer_names(num)) // 'wtl' + ExtDiag(idx)%desc = trim(Model%tracer_names(num)) // ' Large-Scale Wet Deposition' + ExtDiag(idx)%unit = 'kg/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%wetdpl(:,idt) + enddo + + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = trim(Model%tracer_names(num)) // 'wtc' + ExtDiag(idx)%desc = trim(Model%tracer_names(num)) // ' Convective-Scale Wet Deposition' + ExtDiag(idx)%unit = 'kg/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%wetdpc(:,idt) + enddo + endif + enddo + endif + + endif + + num = size(IntDiag(1)%abem, dim=2) + do num = 1, size(IntDiag(1)%abem, dim=2) + idx = idx + 1 + select case (mod(num,3)) + case (0) + ExtDiag(idx)%name = 'bcem' + ExtDiag(idx)%desc = 'Black Carbon' + case (1) + ExtDiag(idx)%name = 'ocem' + ExtDiag(idx)%desc = 'Organic Carbon' + case (2) + ExtDiag(idx)%name = 'so2em' + ExtDiag(idx)%desc = 'SO2' + end select + + if (num > 3) then + ExtDiag(idx)%name = trim(ExtDiag(idx)%name) // 'bb' + ExtDiag(idx)%desc = trim(ExtDiag(idx)%desc) // ' Biomass Burning Emissions' + else + ExtDiag(idx)%name = trim(ExtDiag(idx)%name) // 'an' + ExtDiag(idx)%desc = trim(ExtDiag(idx)%desc) // ' Anthropogenic Emissions' + end if + + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%unit = 'ug/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%abem(:,num) + enddo + end do + + do num = 1, size(IntDiag(1)%aecm, dim=2) + idx = idx + 1 + select case (num) + case(1) + ExtDiag(idx)%name = 'aecmass' + ExtDiag(idx)%desc = 'Aerosol Column Mass Density (PM2.5)' + case(2) + ExtDiag(idx)%name = 'bccmass' + ExtDiag(idx)%desc = 'Black Carbon Column Mass Density' + case(3) + ExtDiag(idx)%name = 'occmass' + ExtDiag(idx)%desc = 'Organic Carbon Column Mass Density' + case(4) + ExtDiag(idx)%name = 'sucmass' + ExtDiag(idx)%desc = 'Sulfate Column Mass Density' + case(5) + ExtDiag(idx)%name = 'ducmass' + ExtDiag(idx)%desc = 'Dust Column Mass Density' + case(6) + ExtDiag(idx)%name = 'sscmass' + ExtDiag(idx)%desc = 'Seasalt Column Mass Density' + end select + + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%unit = 'g/m2' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%aecm(:,num) + enddo + end do + + endif + !--- prognostic variable tendencies (t, u, v, sph, clwmr, o3) !rab idx = idx + 1 !rab ExtDiag(idx)%axes = 3 diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index 1069e3763..51cf35991 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -658,6 +658,10 @@ subroutine GFS_physics_driver & endif endif + if (imp_physics == 99) then + if (Model%cplchm) nvdiff = 3 + end if + ntkev = nvdiff nsteps_per_reset=nint(Model%avg_max_length/dtp) kdtminus1=kdt-1 @@ -1632,6 +1636,17 @@ subroutine GFS_physics_driver & enddo enddo ntiwx = 3 + elseif (imp_physics == 99) then +! Zhao/Carr/Sundqvist + if (Model%cplchm) then + do k=1,levs + do i=1,im + vdftra(i,k,1) = Statein%qgrs(i,k,1) + vdftra(i,k,2) = Statein%qgrs(i,k,ntcw) + vdftra(i,k,3) = Statein%qgrs(i,k,ntoz) + enddo + enddo + endif endif ! if (Model%satmedmf) then @@ -1760,6 +1775,17 @@ subroutine GFS_physics_driver & dqdt(i,k,ntoz) = dvdftra(i,k,7) enddo enddo + + elseif (imp_physics == 99) then + if (Model%cplchm) then + do k=1,levs + do i=1,im + dqdt(i,k,1) = dvdftra(i,k,1) + dqdt(i,k,ntcw) = dvdftra(i,k,2) + dqdt(i,k,ntoz) = dvdftra(i,k,3) + enddo + enddo + endif endif ! if (Model%satmedmf) then @@ -1933,7 +1959,7 @@ subroutine GFS_physics_driver & if (Model%ldiag3d) then do k=1,levs do i=1,im - Diag%dt3dt(i,k,7) = Diag%dt3dt(i,k,7) - dtdt(i,k) + Diag%dt3dt(i,k,7) = Diag%dt3dt(i,k,7) - dtdt(i,k)*dtf enddo enddo endif @@ -3290,14 +3316,14 @@ subroutine GFS_physics_driver & ! enddo ! enddo ! endif - if (Model%ldiag3d) then - do k=1,levs - do i=1,im - Diag%dt3dt(i,k,8) = Diag%dt3dt(i,k,8) + (Stateout%gt0(i,k) -dtdt(i,k) ) * frain -! Diag%dq3dt(i,k,2) = Diag%dq3dt(i,k,2) + (Stateout%gq0(i,k,1)-dqdt(i,k,1)) * frain - enddo - enddo - endif +! if (Model%ldiag3d) then +! do k=1,levs +! do i=1,im +! Diag%dt3dt(i,k,8) = Diag%dt3dt(i,k,8) + (Stateout%gt0(i,k) -dtdt(i,k) ) * frain +!! Diag%dq3dt(i,k,2) = Diag%dq3dt(i,k,2) + (Stateout%gq0(i,k,1)-dqdt(i,k,1)) * frain +! enddo +! enddo +! endif endif endif ! moist convective adjustment over ! @@ -3960,7 +3986,7 @@ subroutine GFS_physics_driver & ! --- ... coupling insertion - if (Model%cplflx) then + if (Model%cplflx .or. Model%cplchm) then do i = 1, im if (t850(i) > 273.16) then Coupling%rain_cpl(i) = Coupling%rain_cpl(i) + Diag%rain(i) @@ -3970,12 +3996,12 @@ subroutine GFS_physics_driver & enddo endif - if (Model%cplchm.and. .not. Model%cplflx) then + if (Model%cplchm) then do i = 1, im - Coupling%rain_cpl(i) = Coupling%rain_cpl(i) + Diag%rain(i) Coupling%rainc_cpl(i) = Coupling%rainc_cpl(i) + Diag%rainc(i) enddo endif + ! --- ... end coupling insertion !!! update surface diagnosis fields at the end of phys package diff --git a/gfsphysics/GFS_layer/GFS_typedefs.F90 b/gfsphysics/GFS_layer/GFS_typedefs.F90 index 20d2b790a..68dcda8d4 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -685,6 +685,9 @@ module GFS_typedefs integer :: nto2 !< tracer index for oxygen integer :: ntwa !< tracer index for water friendly aerosol integer :: ntia !< tracer index for ice friendly aerosol + integer :: ntchm !< number of chemical tracers + integer :: ntchs !< tracer index for first chemical tracer + logical, pointer :: ntdiag(:) => null() !< array to control diagnostics for chemical tracers !--- derived totals for phy_f*d integer :: ntot2d !< total number of variables for phyf2d @@ -1018,10 +1021,22 @@ module GFS_typedefs !--- MP quantities for 3D diagnositics real (kind=kind_phys), pointer :: refl_10cm(:,:) => null() !< instantaneous refl_10cm + !--- Output diagnostics for coupled chemistry + real (kind=kind_phys), pointer :: duem (:,:) => null() !< instantaneous dust emission flux ( kg/m**2/s ) + real (kind=kind_phys), pointer :: ssem (:,:) => null() !< instantaneous sea salt emission flux ( kg/m**2/s ) + real (kind=kind_phys), pointer :: sedim (:,:) => null() !< instantaneous sedimentation ( kg/m**2/s ) + real (kind=kind_phys), pointer :: drydep(:,:) => null() !< instantaneous dry deposition ( kg/m**2/s ) + real (kind=kind_phys), pointer :: wetdpl(:,:) => null() !< instantaneous large-scale wet deposition ( kg/m**2/s ) + real (kind=kind_phys), pointer :: wetdpc(:,:) => null() !< instantaneous convective-scale wet deposition ( kg/m**2/s ) + real (kind=kind_phys), pointer :: abem (:,:) => null() !< instantaneous anthopogenic and biomass burning emissions + !< for black carbon, organic carbon, and sulfur dioxide ( ug/m**2/s ) + real (kind=kind_phys), pointer :: aecm (:,:) => null() !< instantaneous aerosol column mass densities for + !< pm2.5, black carbon, organic carbon, sulfate, dust, sea salt ( g/m**2 ) contains procedure :: create => diag_create procedure :: rad_zero => diag_rad_zero procedure :: phys_zero => diag_phys_zero + procedure :: chem_init => diag_chem_init end type GFS_diag_type #ifdef CCPP @@ -1337,10 +1352,9 @@ subroutine coupling_create (Coupling, IM, Model) Coupling%sfcnsw = clear_val Coupling%sfcdlw = clear_val - if (Model%cplflx .or. Model%do_sppt) then + if (Model%cplflx .or. Model%do_sppt .or. Model%cplchm) then allocate (Coupling%rain_cpl (IM)) allocate (Coupling%snow_cpl (IM)) - Coupling%rain_cpl = clear_val Coupling%snow_cpl = clear_val endif @@ -1490,17 +1504,14 @@ subroutine coupling_create (Coupling, IM, Model) ! -- GSDCHEM coupling options if (Model%cplchm) then !--- outgoing instantaneous quantities - allocate (Coupling%ushfsfci(IM)) - allocate (Coupling%dkt (IM,Model%levs)) - !--- accumulated total and convective rainfall - allocate (Coupling%rain_cpl (IM)) - allocate (Coupling%rainc_cpl (IM)) - - Coupling%rain_cpl = clear_val - Coupling%rainc_cpl = clear_val - - Coupling%ushfsfci = clear_val - Coupling%dkt = clear_val + allocate (Coupling%ushfsfci (IM)) + allocate (Coupling%dkt (IM,Model%levs)) + !--- accumulated convective rainfall + allocate (Coupling%rainc_cpl (IM)) + + Coupling%rainc_cpl = clear_val + Coupling%ushfsfci = clear_val + Coupling%dkt = clear_val endif !--- stochastic physics option @@ -2249,6 +2260,24 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%ntke = get_tracer_index(Model%tracer_names, 'sgs_tke', Model%me, Model%master, Model%debug) Model%ntwa = get_tracer_index(Model%tracer_names, 'liq_aero', Model%me, Model%master, Model%debug) Model%ntia = get_tracer_index(Model%tracer_names, 'ice_aero', Model%me, Model%master, Model%debug) + Model%ntchm = 0 + Model%ntchs = get_tracer_index(Model%tracer_names, 'so2', Model%me, Model%master, Model%debug) + if (Model%ntchs > 0) then + Model%ntchm = get_tracer_index(Model%tracer_names, 'pp10', Model%me, Model%master, Model%debug) + if (Model%ntchm > 0) then + Model%ntchm = Model%ntchm - Model%ntchs + 1 + allocate(Model%ntdiag(Model%ntchm)) + ! -- turn on all tracer diagnostics to .true. by default, except for so2 + Model%ntdiag(1) = .false. + Model%ntdiag(2:) = .true. + ! -- turn off diagnostics for DMS + n = get_tracer_index(Model%tracer_names, 'DMS', Model%me, Model%master, Model%debug) - Model%ntchs + 1 + if (n > 0) Model%ntdiag(n) = .false. + ! -- turn off diagnostics for msa + n = get_tracer_index(Model%tracer_names, 'msa', Model%me, Model%master, Model%debug) - Model%ntchs + 1 + if (n > 0) Model%ntdiag(n) = .false. + endif + endif !--- quantities to be used to derive phy_f*d totals Model%nshoc_2d = nshoc_2d @@ -2844,6 +2873,8 @@ subroutine control_print(Model) print *, ' nto2 : ', Model%nto2 print *, ' ntwa : ', Model%ntwa print *, ' ntia : ', Model%ntia + print *, ' ntchm : ', Model%ntchm + print *, ' ntchs : ', Model%ntchs print *, ' ' print *, 'derived totals for phy_f*d' print *, ' ntot2d : ', Model%ntot2d @@ -3216,15 +3247,21 @@ subroutine diag_create (Diag, IM, Model) ! allocate (Diag%det_mf (IM,Model%levs)) ! allocate (Diag%cldcov (IM,Model%levs)) endif + !--- 3D diagnostics for Thompson MP - allocate (Diag%refl_10cm(IM,Model%levs)) + allocate (Diag%refl_10cm(IM,Model%levs)) + !-- New max hourly diag. - allocate (Diag%refdmax(IM)) - allocate (Diag%refdmax263k(IM)) - allocate (Diag%t02max(IM)) - allocate (Diag%t02min(IM)) - allocate (Diag%rh02max(IM)) - allocate (Diag%rh02min(IM)) + allocate (Diag%refdmax(IM)) + allocate (Diag%refdmax263k(IM)) + allocate (Diag%t02max(IM)) + allocate (Diag%t02min(IM)) + allocate (Diag%rh02max(IM)) + allocate (Diag%rh02min(IM)) + + !--- diagnostics for coupled chemistry + if (Model%cplchm) call Diag%chem_init(IM,Model) + call Diag%rad_zero (Model) ! print *,'in diag_create, call phys_zero' linit = .true. @@ -3379,6 +3416,100 @@ subroutine diag_phys_zero (Diag, Model, linit) endif end subroutine diag_phys_zero +!----------------------- +! GFS_diag%chem_init +!----------------------- + subroutine diag_chem_init(Diag, IM, Model) + + use parse_tracers, only: get_tracer_index, NO_TRACER + + class(GFS_diag_type) :: Diag + integer, intent(in) :: IM + type(GFS_control_type), intent(in) :: Model + + ! -- local variables + integer :: n + + ! -- initialize diagnostic variables depending on + ! -- specific chemical tracers + if (Model%ntchm > 0) then + ! -- retrieve number of dust bins + n = get_number_bins('dust') + if (n > 0) then + allocate (Diag%duem(IM,n)) + Diag%duem = zero + end if + + ! -- retrieve number of sea salt bins + n = get_number_bins('seas') + if (n > 0) then + allocate (Diag%ssem(IM,n)) + Diag%ssem = zero + end if + end if + + ! -- sedimentation and dry/wet deposition diagnostics + if (associated(Model%ntdiag)) then + ! -- get number of tracers with enabled diagnostics + n = count(Model%ntdiag) + + ! -- initialize sedimentation + allocate (Diag%sedim(IM,n)) + Diag%sedim = zero + + ! -- initialize dry deposition + allocate (Diag%drydep(IM,n)) + Diag%drydep = zero + + ! -- initialize large-scale wet deposition + allocate (Diag%wetdpl(IM,n)) + Diag%wetdpl = zero + + ! -- initialize convective-scale wet deposition + allocate (Diag%wetdpc(IM,n)) + Diag%wetdpc = zero + end if + + ! -- initialize anthropogenic and biomass + ! -- burning emission diagnostics for + ! -- (in order): black carbon, + ! -- organic carbon, and sulfur dioxide + allocate (Diag%abem(IM,6)) + Diag%abem = zero + + ! -- initialize column burden diagnostics + ! -- for aerosol species (in order): pm2.5 + ! -- black carbon, organic carbon, sulfate, + ! -- dust, sea salt + allocate (Diag%aecm(IM,6)) + Diag%aecm = zero + + contains + + integer function get_number_bins(tracer_type) + character(len=*), intent(in) :: tracer_type + + logical :: next + integer :: n + character(len=5) :: name + + get_number_bins = 0 + + n = 0 + next = .true. + do while (next) + n = n + 1 + write(name,'(a,i1)') tracer_type, n + 1 + next = get_tracer_index(Model%tracer_names, name, & + Model%me, Model%master, Model%debug) /= NO_TRACER + end do + + get_number_bins = n + + end function get_number_bins + + end subroutine diag_chem_init + #ifdef CCPP ! DH* for testing of CCPP integration !-------------------- diff --git a/gfsphysics/physics/gfdl_cloud_microphys.F90 b/gfsphysics/physics/gfdl_cloud_microphys.F90 index d855fb0a5..609b41ef9 100644 --- a/gfsphysics/physics/gfdl_cloud_microphys.F90 +++ b/gfsphysics/physics/gfdl_cloud_microphys.F90 @@ -840,8 +840,8 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & if (fix_negative) & call neg_adj (ktop, kbot, tz, dp1, qvz, qlz, qrz, qiz, qsz, qgz) - m2_rain (:, :) = 0. - m2_sol (:, :) = 0. + m2_rain (i, :) = 0. + m2_sol (i, :) = 0. do n = 1, ntimes @@ -923,6 +923,10 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & enddo + ! convert units from Pa*kg/kg to kg/m^2/s + m2_rain (i, :) = m2_rain (i, :) * rdt * rgrav + m2_sol (i, :) = m2_sol (i, :) * rdt * rgrav + ! ----------------------------------------------------------------------- ! momentum transportation during sedimentation ! note: dp1 is dry mass; dp0 is the old moist (total) mass diff --git a/module_fcst_grid_comp.F90 b/module_fcst_grid_comp.F90 index f935ba33c..b890de517 100644 --- a/module_fcst_grid_comp.F90 +++ b/module_fcst_grid_comp.F90 @@ -97,14 +97,21 @@ module module_fcst_grid_comp !----- coupled model data ----- integer :: date_init(6) - integer :: numLevels = 0 - integer :: numTracers = 0 + integer :: numLevels = 0 integer :: numSoilLayers = 0 + integer :: numTracers = 0 + integer :: num_diag_sfc_emis_flux = 0 + integer :: num_diag_down_flux = 0 + integer :: num_diag_type_down_flux = 0 + integer :: num_diag_burn_emis_flux = 0 + integer :: num_diag_cmass = 0 ! !----------------------------------------------------------------------- ! public SetServices, fcstGrid - public numLevels, numTracers, numSoilLayers + public numLevels, numSoilLayers, numTracers, & + num_diag_sfc_emis_flux, num_diag_down_flux, & + num_diag_type_down_flux, num_diag_burn_emis_flux, num_diag_cmass ! contains ! @@ -666,8 +673,12 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) !end qulting endif - call get_atmos_model_ungridded_dim(nlev=numLevels, ntracers=numTracers, & - nsoillev=numSoilLayers) + call get_atmos_model_ungridded_dim(nlev=numLevels, nsoillev=numSoilLayers, & + ntracers=numTracers, num_diag_burn_emis_flux=num_diag_burn_emis_flux, & + num_diag_sfc_emis_flux=num_diag_sfc_emis_flux, & + num_diag_down_flux=num_diag_down_flux, & + num_diag_type_down_flux=num_diag_type_down_flux, & + num_diag_cmass=num_diag_cmass) ! !----------------------------------------------------------------------- ! diff --git a/namphysics/NAM_layer/NAM_diagnostics.F90 b/namphysics/NAM_layer/NAM_diagnostics.F90 index 8a4b2bd60..1813bae17 100644 --- a/namphysics/NAM_layer/NAM_diagnostics.F90 +++ b/namphysics/NAM_layer/NAM_diagnostics.F90 @@ -98,7 +98,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop type(GFS_init_type), intent(in) :: Init_parm !--- local variables - integer :: idx, num, nb, nblks, NFXR + integer :: idt, idx, num, nb, nblks, NFXR character(len=2) :: xtra real(kind=kind_phys), parameter :: cn_one = 1._kind_phys real(kind=kind_phys), parameter :: cn_100 = 100._kind_phys @@ -2555,6 +2555,160 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop endif ! print *,'in gfdl_diag_register,af all extdiag, idx=',idx +! -- chemistry diagnostic variables + if (Model%cplchm) then + + if (Model%ntchm > 0) then + + if (associated(IntDiag(1)%duem)) then + do num = 1, size(IntDiag(1)%duem, dim=2) + idx = idx + 1 + ExtDiag(idx)%axes = 2 + write(ExtDiag(idx)%name,'("duem",i3.3)') num + write(ExtDiag(idx)%desc,'("Dust Emission Bin ",i0)') num + ExtDiag(idx)%unit = 'kg/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%duem(:,num) + enddo + enddo + endif + + if (associated(IntDiag(1)%ssem)) then + do num = 1, size(IntDiag(1)%ssem, dim=2) + idx = idx + 1 + ExtDiag(idx)%axes = 2 + write(ExtDiag(idx)%name,'("ssem",i3.3)') num + write(ExtDiag(idx)%desc,'("Seasalt Emission Bin ",i0)') num + ExtDiag(idx)%unit = 'kg/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%ssem(:,num) + enddo + enddo + endif + + if (associated(Model%ntdiag)) then + idt = 0 + do num = Model%ntchs, Model%ntchm + Model%ntchs - 1 + if (Model%ntdiag(num-Model%ntchs+1)) then + idt = idt + 1 + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = trim(Model%tracer_names(num)) // 'sd' + ExtDiag(idx)%desc = trim(Model%tracer_names(num)) // ' Sedimentation' + ExtDiag(idx)%unit = 'kg/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%sedim(:,idt) + enddo + + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = trim(Model%tracer_names(num)) // 'dp' + ExtDiag(idx)%desc = trim(Model%tracer_names(num)) // ' Dry Deposition' + ExtDiag(idx)%unit = 'kg/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%drydep(:,idt) + enddo + + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = trim(Model%tracer_names(num)) // 'wtl' + ExtDiag(idx)%desc = trim(Model%tracer_names(num)) // ' Large-Scale Wet Deposition' + ExtDiag(idx)%unit = 'kg/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%wetdpl(:,idt) + enddo + + idx = idx + 1 + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%name = trim(Model%tracer_names(num)) // 'wtc' + ExtDiag(idx)%desc = trim(Model%tracer_names(num)) // ' Convective-Scale Wet Deposition' + ExtDiag(idx)%unit = 'kg/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%wetdpc(:,idt) + enddo + endif + enddo + endif + + endif + + num = size(IntDiag(1)%abem, dim=2) + do num = 1, size(IntDiag(1)%abem, dim=2) + idx = idx + 1 + select case (mod(num,3)) + case (0) + ExtDiag(idx)%name = 'bcem' + ExtDiag(idx)%desc = 'Black Carbon' + case (1) + ExtDiag(idx)%name = 'ocem' + ExtDiag(idx)%desc = 'Organic Carbon' + case (2) + ExtDiag(idx)%name = 'so2em' + ExtDiag(idx)%desc = 'SO2' + end select + + if (num > 3) then + ExtDiag(idx)%name = trim(ExtDiag(idx)%name) // 'bb' + ExtDiag(idx)%desc = trim(ExtDiag(idx)%desc) // ' Biomass Burning Emissions' + else + ExtDiag(idx)%name = trim(ExtDiag(idx)%name) // 'an' + ExtDiag(idx)%desc = trim(ExtDiag(idx)%desc) // ' Anthropogenic Emissions' + end if + + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%unit = 'ug/m2/s' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%abem(:,num) + enddo + end do + + do num = 1, size(IntDiag(1)%aecm, dim=2) + idx = idx + 1 + select case (num) + case(1) + ExtDiag(idx)%name = 'aecmass' + ExtDiag(idx)%desc = 'Aerosol Column Mass Density (PM2.5)' + case(2) + ExtDiag(idx)%name = 'bccmass' + ExtDiag(idx)%desc = 'Black Carbon Column Mass Density' + case(3) + ExtDiag(idx)%name = 'occmass' + ExtDiag(idx)%desc = 'Organic Carbon Column Mass Density' + case(4) + ExtDiag(idx)%name = 'sucmass' + ExtDiag(idx)%desc = 'Sulfate Column Mass Density' + case(5) + ExtDiag(idx)%name = 'ducmass' + ExtDiag(idx)%desc = 'Dust Column Mass Density' + case(6) + ExtDiag(idx)%name = 'sscmass' + ExtDiag(idx)%desc = 'Seasalt Column Mass Density' + end select + + ExtDiag(idx)%axes = 2 + ExtDiag(idx)%unit = 'g/m2' + ExtDiag(idx)%mod_name = 'gfs_phys' + allocate (ExtDiag(idx)%data(nblks)) + do nb = 1,nblks + ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%aecm(:,num) + enddo + end do + + endif !--- prognostic variable tendencies (t, u, v, sph, clwmr, o3) !rab idx = idx + 1 !rab ExtDiag(idx)%axes = 3 diff --git a/namphysics/NAM_layer/NAM_typedefs.F90 b/namphysics/NAM_layer/NAM_typedefs.F90 index b0efb8c98..6f41455b6 100644 --- a/namphysics/NAM_layer/NAM_typedefs.F90 +++ b/namphysics/NAM_layer/NAM_typedefs.F90 @@ -670,6 +670,9 @@ module GFS_typedefs integer :: nto2 !< tracer index for oxygen integer :: ntwa !< tracer index for water friendly aerosol integer :: ntia !< tracer index for ice friendly aerosol + integer :: ntchm !< number of chemical tracers + integer :: ntchs !< tracer index for first chemical tracer + logical, pointer :: ntdiag(:) => null() !< array to control diagnostics for chemical tracers !--- derived totals for phy_f*d integer :: ntot2d !< total number of variables for phyf2d @@ -1132,10 +1135,22 @@ module GFS_typedefs !--- MP quantities for 3D diagnositics real (kind=kind_phys), pointer :: refl_10cm(:,:) => null() !< instantaneous refl_10cm + !--- Output diagnostics for coupled chemistry + real (kind=kind_phys), pointer :: duem (:,:) => null() !< instantaneous dust emission flux ( kg/m**2/s ) + real (kind=kind_phys), pointer :: ssem (:,:) => null() !< instantaneous sea salt emission flux ( kg/m**2/s ) + real (kind=kind_phys), pointer :: sedim (:,:) => null() !< instantaneous sedimentation ( kg/m**2/s ) + real (kind=kind_phys), pointer :: drydep(:,:) => null() !< instantaneous dry deposition ( kg/m**2/s ) + real (kind=kind_phys), pointer :: wetdpl(:,:) => null() !< instantaneous large-scale wet deposition ( kg/m**2/s ) + real (kind=kind_phys), pointer :: wetdpc(:,:) => null() !< instantaneous convective-scale wet deposition ( kg/m**2/s ) + real (kind=kind_phys), pointer :: abem (:,:) => null() !< instantaneous anthopogenic and biomass burning emissions + !< for black carbon, organic carbon, and sulfur dioxide ( ug/m**2/s ) + real (kind=kind_phys), pointer :: aecm (:,:) => null() !< instantaneous aerosol column mass densities for + !< pm2.5, black carbon, organic carbon, sulfate, dust, sea salt ( g/m**2 ) contains procedure create => diag_create procedure rad_zero => diag_rad_zero procedure phys_zero => diag_phys_zero + procedure chem_init => diag_chem_init end type GFS_diag_type !---------------- @@ -1479,7 +1494,7 @@ subroutine coupling_create (Coupling, IM, Model) Coupling%sfcnsw = clear_val Coupling%sfcdlw = clear_val - if (Model%cplflx .or. Model%do_sppt) then + if (Model%cplflx .or. Model%do_sppt .or. Model%cplchm) then allocate (Coupling%rain_cpl (IM)) allocate (Coupling%snow_cpl (IM)) @@ -1632,17 +1647,14 @@ subroutine coupling_create (Coupling, IM, Model) ! -- GSDCHEM coupling options if (Model%cplchm) then !--- outgoing instantaneous quantities - allocate (Coupling%ushfsfci(IM)) - allocate (Coupling%dkt (IM,Model%levs)) - !--- accumulated total and convective rainfall - allocate (Coupling%rain_cpl (IM)) - allocate (Coupling%rainc_cpl (IM)) - - Coupling%rain_cpl = clear_val - Coupling%rainc_cpl = clear_val - - Coupling%ushfsfci = clear_val - Coupling%dkt = clear_val + allocate (Coupling%ushfsfci (IM)) + allocate (Coupling%dkt (IM,Model%levs)) + !--- accumulated convective rainfall + allocate (Coupling%rainc_cpl (IM)) + + Coupling%rainc_cpl = clear_val + Coupling%ushfsfci = clear_val + Coupling%dkt = clear_val endif !--- stochastic physics option @@ -2423,6 +2435,24 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%ntke = get_tracer_index(Model%tracer_names, 'sgs_tke', Model%me, Model%master, Model%debug) Model%ntwa = get_tracer_index(Model%tracer_names, 'liq_aero', Model%me, Model%master, Model%debug) Model%ntia = get_tracer_index(Model%tracer_names, 'ice_aero', Model%me, Model%master, Model%debug) + Model%ntchm = 0 + Model%ntchs = get_tracer_index(Model%tracer_names, 'so2', Model%me, Model%master, Model%debug) + if (Model%ntchs > 0) then + Model%ntchm = get_tracer_index(Model%tracer_names, 'pp10', Model%me, Model%master, Model%debug) + if (Model%ntchm > 0) then + Model%ntchm = Model%ntchm - Model%ntchs + 1 + allocate(Model%ntdiag(Model%ntchm)) + ! -- turn on all tracer diagnostics to .true. by default, except for so2 + Model%ntdiag(1) = .false. + Model%ntdiag(2:) = .true. + ! -- turn off diagnostics for DMS + n = get_tracer_index(Model%tracer_names, 'DMS', Model%me, Model%master, Model%debug) - Model%ntchs + 1 + if (n > 0) Model%ntdiag(n) = .false. + ! -- turn off diagnostics for msa + n = get_tracer_index(Model%tracer_names, 'msa', Model%me, Model%master, Model%debug) - Model%ntchs + 1 + if (n > 0) Model%ntdiag(n) = .false. + endif + endif !--- quantities to be used to derive phy_f*d totals Model%nshoc_2d = nshoc_2d @@ -3279,6 +3309,8 @@ subroutine control_print(Model) print *, ' nto2 : ', Model%nto2 print *, ' ntwa : ', Model%ntwa print *, ' ntia : ', Model%ntia + print *, ' ntchm : ', Model%ntchm + print *, ' ntchs : ', Model%ntchs print *, ' ' print *, 'derived totals for phy_f*d' print *, ' ntot2d : ', Model%ntot2d @@ -3688,6 +3720,9 @@ subroutine diag_create (Diag, IM, Model) allocate (Diag%refl_10cm(IM,Model%levs)) endif + !--- diagnostics for coupled chemistry + if (Model%cplchm) call Diag%chem_init(IM,Model) + call Diag%rad_zero (Model) ! print *,'in diag_create, call phys_zero' linit = .true. @@ -3903,4 +3938,98 @@ subroutine diag_phys_zero (Diag, Model, linit) endif end subroutine diag_phys_zero +!----------------------- +! GFS_diag%chem_init +!----------------------- + subroutine diag_chem_init(Diag, IM, Model) + + use parse_tracers, only: get_tracer_index, NO_TRACER + + class(GFS_diag_type) :: Diag + integer, intent(in) :: IM + type(GFS_control_type), intent(in) :: Model + + ! -- local variables + integer :: n + + ! -- initialize diagnostic variables depending on + ! -- specific chemical tracers + if (Model%ntchm > 0) then + ! -- retrieve number of dust bins + n = get_number_bins('dust') + if (n > 0) then + allocate (Diag%duem(IM,n)) + Diag%duem = zero + end if + + ! -- retrieve number of sea salt bins + n = get_number_bins('seas') + if (n > 0) then + allocate (Diag%ssem(IM,n)) + Diag%ssem = zero + end if + end if + + ! -- sedimentation and dry/wet deposition diagnostics + if (associated(Model%ntdiag)) then + ! -- get number of tracers with enabled diagnostics + n = count(Model%ntdiag) + + ! -- initialize sedimentation + allocate (Diag%sedim(IM,n)) + Diag%sedim = zero + + ! -- initialize dry deposition + allocate (Diag%drydep(IM,n)) + Diag%drydep = zero + + ! -- initialize large-scale wet deposition + allocate (Diag%wetdpl(IM,n)) + Diag%wetdpl = zero + + ! -- initialize convective-scale wet deposition + allocate (Diag%wetdpc(IM,n)) + Diag%wetdpc = zero + end if + + ! -- initialize anthropogenic and biomass + ! -- burning emission diagnostics for + ! -- (in order): black carbon, + ! -- organic carbon, and sulfur dioxide + allocate (Diag%abem(IM,6)) + Diag%abem = zero + + ! -- initialize column burden diagnostics + ! -- for aerosol species (in order): pm2.5 + ! -- black carbon, organic carbon, sulfate, + ! -- dust, sea salt + allocate (Diag%aecm(IM,6)) + Diag%aecm = zero + + contains + + integer function get_number_bins(tracer_type) + character(len=*), intent(in) :: tracer_type + + logical :: next + integer :: n + character(len=5) :: name + + get_number_bins = 0 + + n = 0 + next = .true. + do while (next) + n = n + 1 + write(name,'(a,i1)') tracer_type, n + 1 + next = get_tracer_index(Model%tracer_names, name, & + Model%me, Model%master, Model%debug) /= NO_TRACER + end do + + get_number_bins = n + + end function get_number_bins + + end subroutine diag_chem_init + end module GFS_typedefs diff --git a/namphysics/physics/gfdl_cloud_microphys.F90 b/namphysics/physics/gfdl_cloud_microphys.F90 index f236dc0d5..57ec6db99 100644 --- a/namphysics/physics/gfdl_cloud_microphys.F90 +++ b/namphysics/physics/gfdl_cloud_microphys.F90 @@ -790,8 +790,8 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & if (fix_negative) & call neg_adj (ktop, kbot, tz, dp1, qvz, qlz, qrz, qiz, qsz, qgz) - m2_rain (:, :) = 0. - m2_sol (:, :) = 0. + m2_rain (i, :) = 0. + m2_sol (i, :) = 0. do n = 1, ntimes @@ -873,6 +873,10 @@ subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, & enddo + ! convert units from Pa*kg/kg to kg/m^2/s + m2_rain (i, :) = m2_rain (i, :) * rdt * rgrav + m2_sol (i, :) = m2_sol (i, :) * rdt * rgrav + ! ----------------------------------------------------------------------- ! momentum transportation during sedimentation ! note: dp1 is dry mass; dp0 is the old moist (total) mass