From 90229c86c9a3c3abd0baefb001d2da979df75ed2 Mon Sep 17 00:00:00 2001 From: wx20hw Date: Thu, 3 Jun 2021 12:46:26 +0000 Subject: [PATCH 01/10] fix the snet bug for Noah LSM --- physics/GFS_radiation_surface.F90 | 2 +- physics/GFS_radiation_surface.meta | 8 +-- physics/gcycle.F90 | 6 ++- physics/radiation_surface.f | 8 ++- physics/sfc_drv.f | 30 ++++++++--- physics/sfc_drv.meta | 81 ++++++++++++++++++++++++++---- 6 files changed, 113 insertions(+), 22 deletions(-) diff --git a/physics/GFS_radiation_surface.F90 b/physics/GFS_radiation_surface.F90 index dd0c56d43..f70413277 100644 --- a/physics/GFS_radiation_surface.F90 +++ b/physics/GFS_radiation_surface.F90 @@ -88,7 +88,7 @@ subroutine GFS_radiation_surface_run ( & alnwf, facsf, facwf, & semis_lnd, semis_ice, snoalb character(len=3) , dimension(:), intent(in) :: lndp_var_list - real(kind=kind_phys), dimension(:), intent(in) :: albdvis_lnd, albdnir_lnd, & + real(kind=kind_phys), dimension(:), intent(inout) :: albdvis_lnd, albdnir_lnd, & albivis_lnd, albinir_lnd real(kind=kind_phys), dimension(:), intent(in) :: albdvis_ice, albdnir_ice, & albivis_ice, albinir_ice diff --git a/physics/GFS_radiation_surface.meta b/physics/GFS_radiation_surface.meta index c38ffe2a3..82faaa865 100644 --- a/physics/GFS_radiation_surface.meta +++ b/physics/GFS_radiation_surface.meta @@ -411,7 +411,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = in + intent = inout optional = F [albdnir_lnd] standard_name = surface_albedo_direct_NIR_over_land @@ -420,7 +420,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = in + intent = inout optional = F [albivis_lnd] standard_name = surface_albedo_diffuse_visible_over_land @@ -429,7 +429,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = in + intent = inout optional = F [albinir_lnd] standard_name = surface_albedo_diffuse_NIR_over_land @@ -438,7 +438,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = in + intent = inout optional = F [albdvis_ice] standard_name = surface_albedo_direct_visible_over_ice diff --git a/physics/gcycle.F90 b/physics/gcycle.F90 index 558a65860..2a33a7e10 100644 --- a/physics/gcycle.F90 +++ b/physics/gcycle.F90 @@ -93,6 +93,7 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, & logical :: lake(nx*ny) + integer :: i_indx(nx*ny), j_indx(nx*ny) character(len=6) :: tile_num_ch real(kind=kind_phys) :: sig1t, dt_warm integer :: npts, nb, ix, jx, ls, ios, ll @@ -120,6 +121,9 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, & end if ! do ix=1,npts + i_indx(ix) = imap(ix) + isc - 1 + j_indx(ix) = jmap(ix) + jsc - 1 + ZORFCS(ix) = zorll (ix) if (slmsk(ix) > 1.9_kind_phys .and. .not. frac_grid) then ZORFCS(ix) = zorli (ix) @@ -190,7 +194,7 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, & nlunit, size(input_nml_file), input_nml_file, & lake, min_lakeice, min_seaice, & ialb, isot, ivegsrc, & - trim(tile_num_ch), imap, jmap) + trim(tile_num_ch), i_indx, j_indx) #ifndef INTERNAL_FILE_NML close (Model%nlunit) #endif diff --git a/physics/radiation_surface.f b/physics/radiation_surface.f index ab7d33e44..b9337636d 100644 --- a/physics/radiation_surface.f +++ b/physics/radiation_surface.f @@ -411,13 +411,15 @@ subroutine setalb & real (kind=kind_phys), dimension(:), intent(in) :: & & slmsk, snowf, zorlf, coszf, tsknf, tairf, hprif, landfrac, & & alvsf, alnsf, alvwf, alnwf, facsf, facwf, fice, tisfc, & - & lsmalbdvis, lsmalbdnir, lsmalbivis, lsmalbinir, & & icealbdvis, icealbdnir, icealbivis, icealbinir, & & sncovr, sncovr_ice, snoalb, albPpert ! sfc-perts, mgehne real (kind=kind_phys), intent(in) :: pertalb ! sfc-perts, mgehne real (kind=kind_phys), intent(in) :: min_seaice real (kind=kind_phys), dimension(:), intent(in) :: & & fracl, fraco, fraci + real (kind=kind_phys), dimension(:),intent(inout) :: & + & lsmalbdvis, lsmalbdnir, lsmalbivis, lsmalbinir + logical, dimension(:), intent(in) :: & & icy @@ -537,6 +539,10 @@ subroutine setalb & alndnd = alnwf(i)*flnd + snoalb(i) * fsno alndvb = ab2bm *flnd + snoalb(i) * fsno alndvd = alvwf(i)*flnd + snoalb(i) * fsno + lsmalbdnir(i) = min(0.99,max(0.01,alndnb)) + lsmalbinir(i) = min(0.99,max(0.01,alndnd)) + lsmalbdvis(i) = min(0.99,max(0.01,alndvb)) + lsmalbivis(i) = min(0.99,max(0.01,alndvd)) else !-- fill in values for land albedo alndnb = 0. diff --git a/physics/sfc_drv.f b/physics/sfc_drv.f index d50a8505e..9ba286095 100644 --- a/physics/sfc_drv.f +++ b/physics/sfc_drv.f @@ -94,10 +94,12 @@ end subroutine lsm_noah_finalize ! call sfc_drv ! ! --- inputs: ! ! ( im, km, ps, t1, q1, soiltyp, vegtype, sigmaf, ! -! sfcemis, dlwflx, dswsfc, snet, delt, tg3, cm, ch, ! +! sfcemis, dlwflx, dswsfc, delt, tg3, cm, ch, ! ! prsl1, prslki, zf, land, wind, slopetyp, ! ! shdmin, shdmax, snoalb, sfalb, flag_iter, flag_guess, ! ! lheatstrg, isot, ivegsrc, ! +! albdvis_lnd, albdnir_lnd, albivis_lnd, albinir_lnd, ! +! adjvisbmd, adjnirbmd, adjvisdfd, adjnirdfd, ! ! --- in/outs: ! ! weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc, ! ! canopy, trans, tsurf, zorl, ! @@ -133,7 +135,6 @@ end subroutine lsm_noah_finalize ! sfcemis - real, sfc lw emissivity ( fraction ) im ! ! dlwflx - real, total sky sfc downward lw flux ( w/m**2 ) im ! ! dswflx - real, total sky sfc downward sw flux ( w/m**2 ) im ! -! snet - real, total sky sfc netsw flx into ground(w/m**2) im ! ! delt - real, time interval (second) 1 ! ! tg3 - real, deep soil temperature (k) im ! ! cm - real, surface exchange coeff for momentum (m/s) im ! @@ -148,6 +149,14 @@ end subroutine lsm_noah_finalize ! shdmax - real, max fractnl cover of green veg (not used) im ! ! snoalb - real, upper bound on max albedo over deep snow im ! ! sfalb - real, mean sfc diffused sw albedo (fractional) im ! +! albdvis_lnd - real, sfc vis direct sw albedo over land im ! +! albdnir_lnd - real, sfc nir direct sw albedo over land im ! +! albivis_lnd - real, sfc vis diffused sw albedo over land im ! +! albinir_lnd - real, sfc nir diffused sw albedo over land im ! +! adjvisbmd - real, t adj sfc uv+vis-beam sw dnward flux im ! +! adjnirbmd - real, t adj sfc nir-beam sw dnward flux im ! +! adjvisdfd - real, t adj sfc uv+vis-diff sw dnward flux im ! +! adjnirdfd - real, t adj sfc nir-diff sw dnward flux im ! ! flag_iter- logical, im ! ! flag_guess-logical, im ! ! lheatstrg- logical, flag for canopy heat storage 1 ! @@ -203,11 +212,13 @@ end subroutine lsm_noah_finalize subroutine lsm_noah_run & & ( im, km, grav, cp, hvap, rd, eps, epsm1, rvrdm1, ps, & ! --- inputs: & t1, q1, soiltyp, vegtype, sigmaf, & - & sfcemis, dlwflx, dswsfc, snet, delt, tg3, cm, ch, & + & sfcemis, dlwflx, dswsfc, delt, tg3, cm, ch, & & prsl1, prslki, zf, land, wind, slopetyp, & & shdmin, shdmax, snoalb, sfalb, flag_iter, flag_guess, & & lheatstrg, isot, ivegsrc, & & bexppert, xlaipert, vegfpert,pertvegf, & ! sfc perts, mgehne + & albdvis_lnd, albdnir_lnd, albivis_lnd, albinir_lnd, & + & adjvisbmd, adjnirbmd, adjvisdfd, adjnirdfd, & ! --- in/outs: & weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc, & & canopy, trans, tsurf, zorl, & @@ -247,10 +258,12 @@ subroutine lsm_noah_run & integer, dimension(:), intent(in) :: soiltyp, vegtype, slopetyp real (kind=kind_phys), dimension(:), intent(in) :: ps, & - & t1, q1, sigmaf, sfcemis, dlwflx, dswsfc, snet, tg3, cm, & + & t1, q1, sigmaf, sfcemis, dlwflx, dswsfc, tg3, cm, & & ch, prsl1, prslki, wind, shdmin, shdmax, & & snoalb, sfalb, zf, & - & bexppert, xlaipert, vegfpert + & bexppert, xlaipert, vegfpert, & + & albdvis_lnd, albdnir_lnd, albivis_lnd, albinir_lnd, & + & adjvisbmd, adjnirbmd, adjvisdfd, adjnirdfd real (kind=kind_phys), intent(in) :: delt @@ -400,7 +413,12 @@ subroutine lsm_noah_run & lwdn = dlwflx(i) !..downward lw flux at sfc in w/m2 swdn = dswsfc(i) !..downward sw flux at sfc in w/m2 - solnet = snet(i) !..net sw rad flx (dn-up) at sfc in w/m2 +!net sw rad flx (dn-up) at sfc in w/m2 + solnet = adjvisbmd(i)*(1-albdvis_lnd(i)) & + & +adjnirbmd(i)*(1-albdnir_lnd(i)) & + & +adjvisdfd(i)*(1-albivis_lnd(i)) & + & +adjnirdfd(i)*(1-albinir_lnd(i)) + sfcems = sfcemis(i) sfcprs = prsl1(i) diff --git a/physics/sfc_drv.meta b/physics/sfc_drv.meta index c68102e7e..7812505c5 100644 --- a/physics/sfc_drv.meta +++ b/physics/sfc_drv.meta @@ -273,15 +273,6 @@ kind = kind_phys intent = in optional = F -[snet] - standard_name = surface_net_downwelling_shortwave_flux - long_name = total sky surface net shortwave flux - units = W m-2 - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in - optional = F [delt] standard_name = time_step_for_dynamics long_name = dynamics time step @@ -482,6 +473,78 @@ kind = kind_phys intent = in optional = F +[albdvis_lnd] + standard_name = surface_albedo_direct_visible_over_land + long_name = direct surface albedo visible band over land + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F +[albdnir_lnd] + standard_name = surface_albedo_direct_NIR_over_land + long_name = direct surface albedo NIR band over land + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F +[albivis_lnd] + standard_name = surface_albedo_diffuse_visible_over_land + long_name = diffuse surface albedo visible band over land + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F +[albinir_lnd] + standard_name = surface_albedo_diffuse_NIR_over_land + long_name = diffuse surface albedo NIR band over land + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F +[adjvisbmd] + standard_name = surface_downwelling_direct_ultraviolet_and_visible_shortwave_flux + long_name = surface downwelling beam ultraviolet plus visible shortwave flux at current time + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F +[adjnirbmd] + standard_name = surface_downwelling_direct_near_infrared_shortwave_flux + long_name = surface downwelling beam near-infrared shortwave flux at current time + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F +[adjvisdfd] + standard_name = surface_downwelling_diffuse_ultraviolet_and_visible_shortwave_flux + long_name = surface downwelling diffuse ultraviolet plus visible shortwave flux at current time + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F +[adjnirdfd] + standard_name = surface_downwelling_diffuse_near_infrared_shortwave_flux + long_name = surface downwelling diffuse near-infrared shortwave flux at current time + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F [weasd] standard_name = water_equivalent_accumulated_snow_depth_over_land long_name = water equiv of acc snow depth over land From d0b52fcad2b7732933851969426e275a4d5a6c49 Mon Sep 17 00:00:00 2001 From: Helin Wei Date: Mon, 7 Jun 2021 23:11:05 -0400 Subject: [PATCH 02/10] include the fix for NoahMP debug test --- physics/module_sf_noahmplsm.f90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index e9e17f399..9fcb7edf8 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -1904,6 +1904,10 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in ! ground snow cover fraction [niu and yang, 2007, jgr] fsno = 0. + if(snowh <= 1.e-6 .or. sneqv <= 1.e-3) then + snowh = 0.0 + sneqv = 0.0 + end if if(snowh.gt.0.) then bdsno = sneqv / snowh fmelt = (bdsno/100.)**parameters%mfsno @@ -7070,7 +7074,7 @@ subroutine snowh2o (parameters,nsnow ,nsoil ,dt ,qsnfro ,qsnsub , & !in end if end if - if(snowh <= 1.e-8 .or. sneqv <= 1.e-6) then + if(snowh <= 1.e-6 .or. sneqv <= 1.e-3) then snowh = 0.0 sneqv = 0.0 end if From 4ccbbeb4b8d244b67493099426bbde87608be9c9 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Sat, 26 Jun 2021 13:13:05 -0600 Subject: [PATCH 03/10] Bugfix in Thompson MP, pass correct timestep to core routine --- physics/mp_thompson.F90 | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 1d235c3e6..0ed8f4a81 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -467,6 +467,9 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & !> - Also, hydrometeor variables are mass or number mixing ratio !> - either kg of species per kg of dry air, or per kg of (dry + vapor). +#if 0 + if (istep==1) then +#endif ! DH* - do this only if istep == 1? Would be ok if it was ! guaranteed that nothing else in the same subcycle group ! was using these arrays, but it is somewhat dangerous. @@ -488,6 +491,9 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if end if ! *DH +#if 0 + endif +#endif !> - Density of air in kg m-3 rho = con_eps*prsl/(con_rd*tgrs*(qv+con_eps)) @@ -565,7 +571,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & if (do_effective_radii) then call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -586,7 +592,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & else call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -607,7 +613,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & else if (do_effective_radii) then call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -627,7 +633,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & first_time_step=first_time_step, errmsg=errmsg, errflg=errflg) else call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -652,6 +658,9 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ! guaranteed that nothing else in the same subcycle group ! was using these arrays, but it is somewhat dangerous. +#if 0 + if(istep==nsteps) then +#endif !> - Convert water vapor mixing ratio back to specific humidity spechum = qv/(1.0_kind_phys+qv) @@ -671,7 +680,9 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if end if ! *DH - +#if 0 + endif +#endif !> - Convert rainfall deltas from mm to m (on physics timestep); add to inout variables ! "rain" in Thompson MP refers to precipitation (total of liquid rainfall+snow+graupel+ice) prcp = prcp + max(0.0, delta_rain_mp/1000.0_kind_phys) From dd04dc3e689854661f460c759bd4efe5991b723e Mon Sep 17 00:00:00 2001 From: Ruiyu Sun Date: Sun, 27 Jun 2021 20:08:59 +0000 Subject: [PATCH 04/10] adding inner loop to Thompson microphysics --- physics/module_mp_thompson.F90 | 97 ++++++++++++++++++++++++++++++---- physics/mp_thompson.F90 | 11 ++-- physics/mp_thompson.meta | 9 ++++ 3 files changed, 103 insertions(+), 14 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 1a038ca72..fa13e574e 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -969,7 +969,7 @@ END SUBROUTINE thompson_init SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & nwfa, nifa, nwfa2d, nifa2d, & tt, th, pii, & - p, w, dz, dt_in, & + p, w, dz, dt_in, dt_inner, & RAINNC, RAINNCV, & SNOWNC, SNOWNCV, & ICENC, ICENCV, & @@ -1029,7 +1029,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & vt_dbz_wt LOGICAL, INTENT(IN) :: first_time_step - REAL, INTENT(IN):: dt_in + REAL, INTENT(IN):: dt_in, dt_inner ! To support subcycling: current step and maximum number of steps INTEGER, INTENT (IN) :: istep, nsteps LOGICAL, INTENT (IN) :: reset_dBZ @@ -1056,6 +1056,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & LOGICAL, OPTIONAL, INTENT(IN) :: diagflag INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref logical :: melti = .false. + INTEGER :: ndt, it ! CCPP error handling character(len=*), optional, intent( out) :: errmsg @@ -1141,7 +1142,25 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & ! j_end = jte ! endif - dt = dt_in +!->rsun minor time interval +! dt = dt_in + RAINNC(:,:) = 0.0 + SNOWNC(:,:) = 0.0 + ICENC(:,:) = 0.0 + GRAUPELNC(:,:) = 0.0 + pcp_ra(:,:) = 0.0 + pcp_sn(:,:) = 0.0 + pcp_gr(:,:) = 0.0 + pcp_ic(:,:) = 0.0 + ndt = max(nint(dt_in/dt_inner),1) + dt = dt_in/ndt + if(dt_in .le. dt_inner) dt= dt_in + if(nsteps > 1) then + write(*,*) 'WARNING: innerloop cant not be used with sybcycle' + ndt = 1 + endif + do it = 1, ndt +!<-rsun qc_max = 0. qr_max = 0. @@ -1260,10 +1279,10 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & rand1, rand2, rand3, & kts, kte, dt, i, j) - pcp_ra(i,j) = pptrain - pcp_sn(i,j) = pptsnow - pcp_gr(i,j) = pptgraul - pcp_ic(i,j) = pptice + pcp_ra(i,j) = pcp_ra(i,j) + pptrain + pcp_sn(i,j) = pcp_sn(i,j) + pptsnow + pcp_gr(i,j) = pcp_gr(i,j) + pptgraul + pcp_ic(i,j) = pcp_ic(i,j) + pptice RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice IF ( PRESENT(snowncv) .AND. PRESENT(snownc) ) THEN @@ -1287,8 +1306,6 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & ENDIF SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12) - - !..Reset lowest model level to initial state aerosols (fake sfc source). !.. Changed 13 May 2013 to fake emissions in which nwfa2d is aerosol !.. number tendency (number per kg per second). @@ -1396,6 +1413,66 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & endif enddo + if (ndt > 1 .and. it == ndt) then + + SR(i,j) = (pcp_sn(i,j) + pcp_gr(i,j) + pcp_ic(i,j))/(RAINNC(i,j)+1.e-12) + RAINNCV(i,j) = RAINNC(i,j) + IF ( PRESENT (snowncv) ) THEN + SNOWNCV(i,j) = SNOWNC(i,j) + ENDIF + IF ( PRESENT (icencv) ) THEN + ICENCV(i,j) = ICENC(i,j) + ENDIF + IF ( PRESENT (graupelncv) ) THEN + GRAUPELNCV(i,j) = GRAUPELNC(i,j) + ENDIF + +!> - Call calc_refl10cm() + + diagflag_presenti: IF ( PRESENT (diagflag) ) THEN + if (diagflag .and. do_radar_ref == 1) then +! + ! Only set melti to true at the output times + if (reset_dBZ) then + melti=.true. + else + melti=.false. + endif +! + if (present(vt_dbz_wt)) then + call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & + t1d, p1d, dBZ, rand1, kts, kte, i, j, & + melti, vt_dbz_wt(i,:,j), & + first_time_step) + else + call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & + t1d, p1d, dBZ, rand1, kts, kte, i, j, & + melti) + end if + do k = kts, kte + refl_10cm(i,k,j) = MAX(-35., dBZ(k)) + enddo + endif + ENDIF diagflag_presenti + + IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN + do k = kts, kte + re_qc1d(k) = re_qc_min + re_qi1d(k) = re_qi_min + re_qs1d(k) = re_qs_min + enddo +!> - Call calc_effectrad() + call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & + re_qc1d, re_qi1d, re_qs1d, kts, kte) + do k = kts, kte + re_cloud(i,k,j) = MAX(re_qc_min, MIN(re_qc1d(k), re_qc_max)) + re_ice(i,k,j) = MAX(re_qi_min, MIN(re_qi1d(k), re_qi_max)) + re_snow(i,k,j) = MAX(re_qs_min, MIN(re_qs1d(k), re_qs_max)) + enddo + ENDIF + + else ! if (ndt > 1 .and. it == ndt) then + ! Diagnostic calculations only for last step ! if Thompson MP is called multiple times last_step_only: IF (istep == nsteps) THEN @@ -1444,6 +1521,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & enddo ENDIF ENDIF last_step_only + endif !if (ndt > 1 .and. it == ndt) then enddo i_loop enddo j_loop @@ -1458,6 +1536,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & ! 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')', & ! 'nr: ', nr_max, '(', imax_nr, ',', jmax_nr, ',', kmax_nr, ')' ! END DEBUG - GT + enddo ! end of nt loop END SUBROUTINE mp_gt_driver !> @} diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 1d235c3e6..57505429f 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -320,7 +320,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & spechum, qc, qr, qi, qs, qg, ni, nr, & is_aerosol_aware, nc, nwfa, nifa, & nwfa2d, nifa2d, & - tgrs, prsl, phii, omega, & + tgrs, prsl, phii, omega, dt_inner, & dtp, first_time_step, istep, nsteps, & prcp, rain, graupel, ice, snow, sr, & refl_10cm, reset_dBZ, do_radar_ref, & @@ -374,6 +374,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ! Radar reflectivity real(kind_phys), intent( out) :: refl_10cm(:,:) logical, optional, intent(in ) :: do_radar_ref + real, intent(in ) :: dt_inner ! Cloud effective radii real(kind_phys), optional, intent( out) :: re_cloud(:,:) real(kind_phys), optional, intent( out) :: re_ice(:,:) @@ -565,7 +566,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & if (do_effective_radii) then call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, dt_inner=dt_inner, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -586,7 +587,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & else call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, dt_inner=dt_inner, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -607,7 +608,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & else if (do_effective_radii) then call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, dt_inner=dt_inner, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -627,7 +628,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & first_time_step=first_time_step, errmsg=errmsg, errflg=errflg) else call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, dt_inner=dt_inner, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta index 573bab6c8..3a3977f88 100644 --- a/physics/mp_thompson.meta +++ b/physics/mp_thompson.meta @@ -545,6 +545,15 @@ kind = kind_phys intent = in optional = F +[dt_inner] + standard_name = time_step_for_inner_loop + long_name = time step for inner loop + units = s + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [first_time_step] standard_name = flag_for_first_time_step long_name = flag for first time step for time integration loop (cold/warmstart) From 1ccfb69bd8fc55e044875629506114cd2e34538c Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 28 Jun 2021 07:10:01 -0600 Subject: [PATCH 05/10] Remove test code from physics/mp_thompson.F90 --- physics/mp_thompson.F90 | 21 +++++---------------- 1 file changed, 5 insertions(+), 16 deletions(-) diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 0ed8f4a81..f4720eb72 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -467,9 +467,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & !> - Also, hydrometeor variables are mass or number mixing ratio !> - either kg of species per kg of dry air, or per kg of (dry + vapor). -#if 0 - if (istep==1) then -#endif ! DH* - do this only if istep == 1? Would be ok if it was ! guaranteed that nothing else in the same subcycle group ! was using these arrays, but it is somewhat dangerous. @@ -491,9 +488,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if end if ! *DH -#if 0 - endif -#endif !> - Density of air in kg m-3 rho = con_eps*prsl/(con_rd*tgrs*(qv+con_eps)) @@ -571,7 +565,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & if (do_effective_radii) then call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -592,7 +586,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & else call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -613,7 +607,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & else if (do_effective_radii) then call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -633,7 +627,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & first_time_step=first_time_step, errmsg=errmsg, errflg=errflg) else call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -658,9 +652,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ! guaranteed that nothing else in the same subcycle group ! was using these arrays, but it is somewhat dangerous. -#if 0 - if(istep==nsteps) then -#endif !> - Convert water vapor mixing ratio back to specific humidity spechum = qv/(1.0_kind_phys+qv) @@ -680,9 +671,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if end if ! *DH -#if 0 - endif -#endif + !> - Convert rainfall deltas from mm to m (on physics timestep); add to inout variables ! "rain" in Thompson MP refers to precipitation (total of liquid rainfall+snow+graupel+ice) prcp = prcp + max(0.0, delta_rain_mp/1000.0_kind_phys) From 7178cb279f8b6975ade28bb248cc7d3108ca79fc Mon Sep 17 00:00:00 2001 From: Ruiyu Sun Date: Tue, 29 Jun 2021 09:52:31 +0000 Subject: [PATCH 06/10] change the logic --- physics/module_mp_thompson.F90 | 83 ++++++++-------------------------- 1 file changed, 19 insertions(+), 64 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index fa13e574e..6ae381195 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1142,7 +1142,6 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & ! j_end = jte ! endif -!->rsun minor time interval ! dt = dt_in RAINNC(:,:) = 0.0 SNOWNC(:,:) = 0.0 @@ -1159,8 +1158,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & write(*,*) 'WARNING: innerloop cant not be used with sybcycle' ndt = 1 endif + do it = 1, ndt -!<-rsun qc_max = 0. qr_max = 0. @@ -1413,73 +1412,30 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & endif enddo - if (ndt > 1 .and. it == ndt) then - - SR(i,j) = (pcp_sn(i,j) + pcp_gr(i,j) + pcp_ic(i,j))/(RAINNC(i,j)+1.e-12) - RAINNCV(i,j) = RAINNC(i,j) - IF ( PRESENT (snowncv) ) THEN - SNOWNCV(i,j) = SNOWNC(i,j) - ENDIF - IF ( PRESENT (icencv) ) THEN - ICENCV(i,j) = ICENC(i,j) - ENDIF - IF ( PRESENT (graupelncv) ) THEN - GRAUPELNCV(i,j) = GRAUPELNC(i,j) - ENDIF - -!> - Call calc_refl10cm() + if (ndt>1 .and. it==ndt) then - diagflag_presenti: IF ( PRESENT (diagflag) ) THEN - if (diagflag .and. do_radar_ref == 1) then -! - ! Only set melti to true at the output times - if (reset_dBZ) then - melti=.true. - else - melti=.false. - endif -! - if (present(vt_dbz_wt)) then - call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & - t1d, p1d, dBZ, rand1, kts, kte, i, j, & - melti, vt_dbz_wt(i,:,j), & - first_time_step) - else - call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & - t1d, p1d, dBZ, rand1, kts, kte, i, j, & - melti) - end if - do k = kts, kte - refl_10cm(i,k,j) = MAX(-35., dBZ(k)) - enddo - endif - ENDIF diagflag_presenti - - IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN - do k = kts, kte - re_qc1d(k) = re_qc_min - re_qi1d(k) = re_qi_min - re_qs1d(k) = re_qs_min - enddo -!> - Call calc_effectrad() - call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & - re_qc1d, re_qi1d, re_qs1d, kts, kte) - do k = kts, kte - re_cloud(i,k,j) = MAX(re_qc_min, MIN(re_qc1d(k), re_qc_max)) - re_ice(i,k,j) = MAX(re_qi_min, MIN(re_qi1d(k), re_qi_max)) - re_snow(i,k,j) = MAX(re_qs_min, MIN(re_qs1d(k), re_qs_max)) - enddo + SR(i,j) = (pcp_sn(i,j) + pcp_gr(i,j) + pcp_ic(i,j))/(RAINNC(i,j)+1.e-12) + RAINNCV(i,j) = RAINNC(i,j) + IF ( PRESENT (snowncv) ) THEN + SNOWNCV(i,j) = SNOWNC(i,j) ENDIF - - else ! if (ndt > 1 .and. it == ndt) then + IF ( PRESENT (icencv) ) THEN + ICENCV(i,j) = ICENC(i,j) + ENDIF + IF ( PRESENT (graupelncv) ) THEN + GRAUPELNCV(i,j) = GRAUPELNC(i,j) + ENDIF + endif ! Diagnostic calculations only for last step ! if Thompson MP is called multiple times - last_step_only: IF (istep == nsteps) THEN + last_step_only:if ((ndt>1 .and. it==ndt) & + .or. (nsteps>1 .and. istep==nsteps) & + .or. (nsteps==1 .and. ndt==1)) then !> - Call calc_refl10cm() - diagflag_present: IF ( PRESENT (diagflag) ) THEN + diagflag_presenti: IF ( PRESENT (diagflag) ) THEN if (diagflag .and. do_radar_ref == 1) then ! ! Only set melti to true at the output times @@ -1503,7 +1459,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & refl_10cm(i,k,j) = MAX(-35., dBZ(k)) enddo endif - ENDIF diagflag_present + ENDIF diagflag_presenti IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN do k = kts, kte @@ -1520,9 +1476,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & re_snow(i,k,j) = MAX(re_qs_min, MIN(re_qs1d(k), re_qs_max)) enddo ENDIF - ENDIF last_step_only - endif !if (ndt > 1 .and. it == ndt) then + endif last_step_only enddo i_loop enddo j_loop From f7450c2bc6194f7f6f2671bbcda99e907ef002bf Mon Sep 17 00:00:00 2001 From: Ruiyu Sun Date: Thu, 1 Jul 2021 04:39:07 +0000 Subject: [PATCH 07/10] remove extra i --- physics/module_mp_thompson.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 6ae381195..92d79a0c1 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1435,7 +1435,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & !> - Call calc_refl10cm() - diagflag_presenti: IF ( PRESENT (diagflag) ) THEN + diagflag_present: IF ( PRESENT (diagflag) ) THEN if (diagflag .and. do_radar_ref == 1) then ! ! Only set melti to true at the output times @@ -1459,7 +1459,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & refl_10cm(i,k,j) = MAX(-35., dBZ(k)) enddo endif - ENDIF diagflag_presenti + ENDIF diagflag_present IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN do k = kts, kte From 581f5a6b7675e142a1c7f8f8d1d22e8a6164cad1 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Sat, 26 Jun 2021 13:13:05 -0600 Subject: [PATCH 08/10] Bugfix in Thompson MP, pass correct timestep to core routine --- physics/mp_thompson.F90 | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 57505429f..1eb2506ae 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -468,6 +468,9 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & !> - Also, hydrometeor variables are mass or number mixing ratio !> - either kg of species per kg of dry air, or per kg of (dry + vapor). +#if 0 + if (istep==1) then +#endif ! DH* - do this only if istep == 1? Would be ok if it was ! guaranteed that nothing else in the same subcycle group ! was using these arrays, but it is somewhat dangerous. @@ -489,6 +492,9 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if end if ! *DH +#if 0 + endif +#endif !> - Density of air in kg m-3 rho = con_eps*prsl/(con_rd*tgrs*(qv+con_eps)) @@ -566,7 +572,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & if (do_effective_radii) then call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, dt_inner=dt_inner, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -587,7 +593,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & else call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, dt_inner=dt_inner, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -608,7 +614,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & else if (do_effective_radii) then call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, dt_inner=dt_inner, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -628,7 +634,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & first_time_step=first_time_step, errmsg=errmsg, errflg=errflg) else call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & - tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, dt_inner=dt_inner, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & icenc=ice_mp, icencv=delta_ice_mp, & @@ -653,6 +659,9 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ! guaranteed that nothing else in the same subcycle group ! was using these arrays, but it is somewhat dangerous. +#if 0 + if(istep==nsteps) then +#endif !> - Convert water vapor mixing ratio back to specific humidity spechum = qv/(1.0_kind_phys+qv) @@ -672,7 +681,9 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if end if ! *DH - +#if 0 + endif +#endif !> - Convert rainfall deltas from mm to m (on physics timestep); add to inout variables ! "rain" in Thompson MP refers to precipitation (total of liquid rainfall+snow+graupel+ice) prcp = prcp + max(0.0, delta_rain_mp/1000.0_kind_phys) From b498507041682c5cbe150f186450db4402d2b196 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 28 Jun 2021 07:10:01 -0600 Subject: [PATCH 09/10] Remove test code from physics/mp_thompson.F90 --- physics/mp_thompson.F90 | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 1eb2506ae..2115615d5 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -468,9 +468,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & !> - Also, hydrometeor variables are mass or number mixing ratio !> - either kg of species per kg of dry air, or per kg of (dry + vapor). -#if 0 - if (istep==1) then -#endif ! DH* - do this only if istep == 1? Would be ok if it was ! guaranteed that nothing else in the same subcycle group ! was using these arrays, but it is somewhat dangerous. @@ -492,9 +489,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if end if ! *DH -#if 0 - endif -#endif !> - Density of air in kg m-3 rho = con_eps*prsl/(con_rd*tgrs*(qv+con_eps)) @@ -659,9 +653,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ! guaranteed that nothing else in the same subcycle group ! was using these arrays, but it is somewhat dangerous. -#if 0 - if(istep==nsteps) then -#endif !> - Convert water vapor mixing ratio back to specific humidity spechum = qv/(1.0_kind_phys+qv) @@ -681,9 +672,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if end if ! *DH -#if 0 - endif -#endif + !> - Convert rainfall deltas from mm to m (on physics timestep); add to inout variables ! "rain" in Thompson MP refers to precipitation (total of liquid rainfall+snow+graupel+ice) prcp = prcp + max(0.0, delta_rain_mp/1000.0_kind_phys) From cdbe3749f365ca178dc0ee533e66f693eef5835a Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Thu, 22 Jul 2021 09:00:54 -0600 Subject: [PATCH 10/10] Remove unneeded block of code from physics/gcycle.F90 --- physics/gcycle.F90 | 7 ------- 1 file changed, 7 deletions(-) diff --git a/physics/gcycle.F90 b/physics/gcycle.F90 index 6431a6954..718b375af 100644 --- a/physics/gcycle.F90 +++ b/physics/gcycle.F90 @@ -172,13 +172,6 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, min_ice(ix) = min_seaice endif - ZORFCS(ix) = zorll (ix) - if (slmsk(ix) > 1.9_kind_phys .and. .not. frac_grid) then - ZORFCS(ix) = zorli (ix) - elseif (slmsk(ix) < 0.1_kind_phys .and. .not. frac_grid) then - ZORFCS(ix) = zorlo (ix) - endif - IF (slmsk(ix) > 1.99_kind_phys) THEN AISFCS(ix) = 1.0_kind_phys ELSE