From c58e27a8641e869e4218bf8626dc6304a32fc31f Mon Sep 17 00:00:00 2001 From: "Shan.Sun" Date: Wed, 8 Jan 2020 21:55:33 +0000 Subject: [PATCH 01/15] Bugfix in the cap at lake points; fice is wrt water area for sea ice and lake ice --- atmos_model.F90 | 11 ++--- gfsphysics/GFS_layer/GFS_physics_driver.F90 | 50 +++++++++------------ io/FV3GFS_io.F90 | 47 +++++++++++-------- 3 files changed, 52 insertions(+), 56 deletions(-) diff --git a/atmos_model.F90 b/atmos_model.F90 index 23e30e76c..af24b1511 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -1697,18 +1697,15 @@ subroutine assign_importdata(rc) ix = Atm_block%ixp(i,j) IPD_Data(nb)%Coupling%ficein_cpl(ix) = zero if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then - if (datar8(i,j) >= IPD_control%min_seaice*IPD_Data(nb)%Sfcprop%oceanfrac(ix)) then + datar8(i,j) = datar8(i,j)/IPD_Data(nb)%Sfcprop%oceanfrac(ix) !LHS: ice frac wrt non-land portion (not whole cell) + if (datar8(i,j) >= IPD_control%min_seaice) then IPD_Data(nb)%Coupling%ficein_cpl(ix) = datar8(i,j) -! if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) == one) IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points - IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points + if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) == one) IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points IPD_Data(nb)%Coupling%slimskin_cpl(ix) = 4. else if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) == one) IPD_Data(nb)%Sfcprop%slmsk(ix) = zero IPD_Data(nb)%Coupling%slimskin_cpl(ix) = zero endif - else - IPD_Data(nb)%Sfcprop%slmsk(ix) = one - IPD_Data(nb)%Coupling%slimskin_cpl(ix) = one endif enddo enddo @@ -1878,7 +1875,7 @@ subroutine assign_importdata(rc) ix = Atm_block%ixp(i,j) if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then !if it is ocean or ice get surface temperature from mediator - if(IPD_Data(nb)%Coupling%ficein_cpl(ix) >= IPD_control%min_seaice*IPD_Data(nb)%Sfcprop%oceanfrac(ix)) then + if(IPD_Data(nb)%Coupling%ficein_cpl(ix) >= IPD_control%min_seaice) then IPD_Data(nb)%Sfcprop%tisfc(ix) = IPD_Data(nb)%Coupling%tisfcin_cpl(ix) IPD_Data(nb)%Sfcprop%fice(ix) = IPD_Data(nb)%Coupling%ficein_cpl(ix) IPD_Data(nb)%Sfcprop%hice(ix) = IPD_Data(nb)%Coupling%hicein_cpl(ix) diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index 5b67f7faa..cf8722765 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -1055,7 +1055,7 @@ subroutine GFS_physics_driver & endif ! --- ... xw: transfer ice thickness & concentration from global to local variables zice(i) = Sfcprop%hice(i) - fice(i) = Sfcprop%fice(i) ! ice fraction of lake/ocean wrt whole cell + fice(i) = min(one,Sfcprop%fice(i)) ! sea/lake ice fraction wrt water portion, not whole cell tice(i) = Sfcprop%tisfc(i) ! !GFDL work1(i) = (log(coslat(i) / (nlons(i)*latr)) - dxmin) * dxinv @@ -1098,22 +1098,20 @@ subroutine GFS_physics_driver & endif ! DH* In CCPP, this is in GFS_surface_composites_pre - if (Model%frac_grid) then ! here Sfcprop%fice is fraction of the whole grid that is ice + if (Model%frac_grid) then do i = 1, IM frland(i) = Sfcprop%landfrac(i) if (frland(i) > zero) dry(i) = .true. - tem = one - frland(i) - if (tem > zero) then + if (frland(i) < one) then if (flag_cice(i)) then - if (fice(i) >= Model%min_seaice*tem) then + if (fice(i) >= Model%min_seaice) then icy(i) = .true. else fice(i) = zero endif else - if (fice(i) >= Model%min_lakeice*tem) then + if (fice(i) >= Model%min_lakeice) then icy(i) = .true. - fice(i) = fice(i)/tem ! fice is fraction of ocean/lake else fice(i) = zero endif @@ -1122,14 +1120,7 @@ subroutine GFS_physics_driver & else fice(i) = zero endif - ! ocean/lake area that is not frozen - tem1 = max(zero, tem - Sfcprop%fice(i)) - - if (tem1 > zero) then - wet(i) = .true. ! there is some open water! -! if (icy(i)) Sfcprop%tsfco(i) = max(Sfcprop%tsfco(i), tgice) - if (icy(i)) Sfcprop%tsfco(i) = max(Sfcprop%tisfc(i), tgice) - endif + if (frland(i)+fice(i)*(one-frland(i)) < one) wet(i)=.true. !there is some open ocean/lake water! enddo else do i = 1, IM @@ -1198,7 +1189,7 @@ subroutine GFS_physics_driver & enddo endif do i=1,im - if(wet(i)) then ! Water + if(wet(i) .or. icy(i)) then ! open water or ice zorl3(i,3) = Sfcprop%zorlo(i) tsfc3(i,3) = Sfcprop%tsfco(i) tsurf3(i,3) = Sfcprop%tsfco(i) @@ -1482,15 +1473,15 @@ subroutine GFS_physics_driver & if (Model%frac_grid) then do i=1,im - tem = one - Sfcprop%fice(i) - frland(i) + tem = (one - frland(i)) * fice(i) ! tem = ice fraction wrt whole cell if (flag_cice(i)) then - adjsfculw(i) = adjsfculw3(i,1) * frland(i) & - + Coupling%ulwsfcin_cpl(i) * Sfcprop%fice(i) & - + adjsfculw3(i,3) * tem + adjsfculw(i) = adjsfculw3(i,1) * frland(i) & + + Coupling%ulwsfcin_cpl(i) * tem & + + adjsfculw3(i,3) * (one - frland(i) - tem) else adjsfculw(i) = adjsfculw3(i,1) * frland(i) & - + adjsfculw3(i,2) * Sfcprop%fice(i) & - + adjsfculw3(i,3) * tem + + adjsfculw3(i,2) * tem & + + adjsfculw3(i,3) * (one - frland(i) - tem) endif enddo else @@ -1498,16 +1489,16 @@ subroutine GFS_physics_driver & if (dry(i)) then ! all land adjsfculw(i) = adjsfculw3(i,1) elseif (icy(i)) then ! ice (and water) - tem = one - Sfcprop%fice(i) + tem = one - fice(i) if (flag_cice(i)) then if (wet(i) .and. adjsfculw3(i,3) /= huge) then - adjsfculw(i) = Coupling%ulwsfcin_cpl(i)*Sfcprop%fice(i) + adjsfculw3(i,3)*tem + adjsfculw(i) = Coupling%ulwsfcin_cpl(i)*fice(i) + adjsfculw3(i,3)*tem else adjsfculw(i) = Coupling%ulwsfcin_cpl(i) endif else if (wet(i) .and. adjsfculw3(i,3) /= huge) then - adjsfculw(i) = adjsfculw3(i,2)*Sfcprop%fice(i) + adjsfculw3(i,3)*tem + adjsfculw(i) = adjsfculw3(i,2)*fice(i) + adjsfculw3(i,3)*tem else adjsfculw(i) = adjsfculw3(i,2) endif @@ -1934,9 +1925,9 @@ subroutine GFS_physics_driver & do i=1, im ! ! Three-way composites (fields from sfc_diff) - txl = Sfcprop%landfrac(i) - txi = Sfcprop%fice(i) ! here Sfcprop%fice is grid fraction that is ice - txo = one - txl - txi + txl = frland(i) + txi = fice(i)*(one - frland(i)) ! txi = ice fraction wrt whole cell + txo = max(zero, one - txl - txi) Sfcprop%zorl(i) = txl*zorl3(i,1) + txi*zorl3(i,2) + txo*zorl3(i,3) cd(i) = txl*cd3(i,1) + txi*cd3(i,2) + txo*cd3(i,3) cdq(i) = txl*cdq3(i,1) + txi*cdq3(i,2) + txo*cdq3(i,3) @@ -1994,8 +1985,7 @@ subroutine GFS_physics_driver & if (.not. flag_cice(i)) then if (islmsk(i) == 2) then ! return updated lake ice thickness & concentration to global array Sfcprop%hice(i) = zice(i) -! Sfcprop%fice(i) = fice(i) * Sfcprop%lakefrac(i) ! fice is fraction of lake area that is frozen - Sfcprop%fice(i) = fice(i) * (one-Sfcprop%landfrac(i)) ! fice is fraction of wet area that is frozen + Sfcprop%fice(i) = fice(i) Sfcprop%tisfc(i) = tice(i) else ! this would be over open ocean or land (no ice fraction) Sfcprop%hice(i) = zero diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index 4b2d1e426..2e19b26b9 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -29,6 +29,7 @@ module FV3GFS_io_mod use diag_data_mod, only: output_fields, max_output_fields use diag_util_mod, only: find_input_field use constants_mod, only: grav, rdgas + use physcons, only: con_tice !saltwater freezing temp (K) ! !--- GFS physics modules !#ifndef CCPP @@ -104,6 +105,7 @@ module FV3GFS_io_mod real, parameter :: missing_value = 9.99e20 real, parameter:: stndrd_atmos_ps = 101325. real, parameter:: stndrd_atmos_lapse = 0.0065 + real, parameter:: drythresh = 1.e-4 !--- miscellaneous other variables logical :: use_wrtgridcomp_output = .FALSE. @@ -1083,7 +1085,8 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) do nb = 1, Atm_block%nblks do ix = 1, Atm_block%blksz(nb) Sfcprop(nb)%sncovr(ix) = 0.0 - if (Sfcprop(nb)%slmsk(ix) > 0.001) then + !if (Sfcprop(nb)%slmsk(ix) > 0.001) then + if (Sfcprop(nb)%landfrac(ix) >= drythresh) then vegtyp = Sfcprop(nb)%vtype(ix) if (vegtyp == 0) vegtyp = 7 rsnow = 0.001*Sfcprop(nb)%weasd(ix)/snupx(vegtyp) @@ -1098,16 +1101,31 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) endif !#endif + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + if (Sfcprop(nb)%lakefrac(ix) > 0.0) then + Sfcprop(nb)%oceanfrac(ix) = 0.0 ! lake & ocean don't coexist in a cell + if (Sfcprop(nb)%fice(ix) < Model%min_lakeice) Sfcprop(nb)%fice(ix) = 0. + else + Sfcprop(nb)%oceanfrac(ix) = 1.0 - Sfcprop(nb)%landfrac(ix) !LHS:ocean frac [0:1] + if (Sfcprop(nb)%fice(ix) < Model%min_seaice) Sfcprop(nb)%fice(ix) = 0. + endif + Sfcprop(nb)%slmsk(ix) = ceiling(Sfcprop(nb)%landfrac(ix)) + if (Sfcprop(nb)%fice(ix) > 0. .and. Sfcprop(nb)%landfrac(ix)==0.) Sfcprop(nb)%slmsk(ix) = 2 ! land dominates over ice if co-exist + enddo + enddo + if(Model%frac_grid) then ! 3-way composite do nb = 1, Atm_block%nblks do ix = 1, Atm_block%blksz(nb) - tem = 1.0 - Sfcprop(nb)%landfrac(ix) - Sfcprop(nb)%fice(ix) + Sfcprop(nb)%tsfco(ix) = max(con_tice, Sfcprop(nb)%tsfco(ix)) + tem = (1.-Sfcprop(nb)%landfrac(ix)) * Sfcprop(nb)%fice(ix) ! tem = ice fraction wrt whole cell Sfcprop(nb)%zorl(ix) = Sfcprop(nb)%zorll(ix) * Sfcprop(nb)%landfrac(ix) & - + Sfcprop(nb)%zorll(ix) * Sfcprop(nb)%fice(ix) & !zorl ice = zorl land - + Sfcprop(nb)%zorlo(ix) * tem + + Sfcprop(nb)%zorll(ix) * tem & !zorl ice = zorl land + + Sfcprop(nb)%zorlo(ix) * (1.-Sfcprop(nb)%landfrac(ix)-tem) Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tsfcl(ix) * Sfcprop(nb)%landfrac(ix) & - + Sfcprop(nb)%tisfc(ix) * Sfcprop(nb)%fice(ix) & - + Sfcprop(nb)%tsfco(ix) * tem + + Sfcprop(nb)%tisfc(ix) * tem & + + Sfcprop(nb)%tsfco(ix) * (1.-Sfcprop(nb)%landfrac(ix)-tem) enddo enddo else ! in this case ice fracion is fraction of water fraction @@ -1118,6 +1136,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) Sfcprop(nb)%zorll(ix) = Sfcprop(nb)%zorlo(ix) Sfcprop(nb)%zorl(ix) = Sfcprop(nb)%zorlo(ix) Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tsfco(ix) + Sfcprop(nb)%tsfco(ix) = max(con_tice, Sfcprop(nb)%tsfco(ix)) if (Sfcprop(nb)%slmsk(ix) > 1.9) then Sfcprop(nb)%landfrac(ix) = 0.0 else @@ -1127,17 +1146,6 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) enddo endif ! if (Model%frac_grid) - do nb = 1, Atm_block%nblks - do ix = 1, Atm_block%blksz(nb) - if (Sfcprop(nb)%lakefrac(ix) > 0.0) then - Sfcprop(nb)%oceanfrac(ix) = 0.0 ! lake & ocean don't coexist in a cell - else - Sfcprop(nb)%oceanfrac(ix) = 1.0 - Sfcprop(nb)%landfrac(ix) !LHS:ocean frac [0:1] - endif - - enddo - enddo - if (Model%lsm == Model%lsm_noahmp) then if (nint(sfc_var2(1,1,nvar_s2m+19)) == -66666) then if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver:: - Cold start Noah MP ') @@ -1182,7 +1190,8 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) Sfcprop(nb)%smoiseq(ix, 1:4) = missing_value Sfcprop(nb)%zsnsoxy(ix, -2:4) = missing_value - if (Sfcprop(nb)%slmsk(ix) > 0.01) then + ! if (Sfcprop(nb)%slmsk(ix) > 0.01) then + if (Sfcprop(nb)%landfrac(ix) >= drythresh) then Sfcprop(nb)%tvxy(ix) = Sfcprop(nb)%tsfcl(ix) Sfcprop(nb)%tgxy(ix) = Sfcprop(nb)%tsfcl(ix) @@ -1385,7 +1394,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) smc = smc - dx if ( abs (dx) < 1.e-6) exit enddo ! iteration - Sfcprop(nb)%smoiseq(ix,ns) = min(max(smc,1.e-4),smcmax*0.99) + Sfcprop(nb)%smoiseq(ix,ns) = min(max(smc,drythresh),smcmax*0.99) enddo ! ddz soil layer else ! bexp <= 0.0 Sfcprop(nb)%smoiseq(ix,1:4) = smcmax From e15f30441ca0112864ba3f52ac3fa1986be41e5a Mon Sep 17 00:00:00 2001 From: "Shan.Sun" Date: Thu, 9 Jan 2020 05:51:17 +0000 Subject: [PATCH 02/15] assign tsfco/zorlo (over water) to tsfcl/zorll (over land) initially, if latter don't exist in sfc_data. --- io/FV3GFS_io.F90 | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index 2e19b26b9..87d23e624 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -1099,6 +1099,29 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) enddo enddo endif + + if(Model%cplflx) then + if (nint(sfc_var2(1,1,33)) == -9999) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing tsfcl') + !--- compute tsfcl from existing variables + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + Sfcprop(nb)%tsfcl(ix) = Sfcprop(nb)%tsfco(ix) + enddo + enddo + endif + + if (nint(sfc_var2(1,1,34)) == -9999) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorll') + !--- compute zorll from existing variables + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + Sfcprop(nb)%zorll(ix) = Sfcprop(nb)%zorlo(ix) + enddo + enddo + endif + endif + !#endif do nb = 1, Atm_block%nblks From 1f61c0b9784cb5b523bdc2eda9134c462e4f66b6 Mon Sep 17 00:00:00 2001 From: "Shan.Sun" Date: Thu, 9 Jan 2020 16:52:55 +0000 Subject: [PATCH 03/15] Limit ice frac to be <=1. --- gfsphysics/GFS_layer/GFS_physics_driver.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index cf8722765..f40f3162d 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -1055,7 +1055,8 @@ subroutine GFS_physics_driver & endif ! --- ... xw: transfer ice thickness & concentration from global to local variables zice(i) = Sfcprop%hice(i) - fice(i) = min(one,Sfcprop%fice(i)) ! sea/lake ice fraction wrt water portion, not whole cell + Sfcprop%fice(i) = min(one,Sfcprop%fice(i)) ! sea/lake ice fraction wrt water portion, not whole cell + fice(i) = Sfcprop%fice(i) tice(i) = Sfcprop%tisfc(i) ! !GFDL work1(i) = (log(coslat(i) / (nlons(i)*latr)) - dxmin) * dxinv @@ -1985,7 +1986,7 @@ subroutine GFS_physics_driver & if (.not. flag_cice(i)) then if (islmsk(i) == 2) then ! return updated lake ice thickness & concentration to global array Sfcprop%hice(i) = zice(i) - Sfcprop%fice(i) = fice(i) + Sfcprop%fice(i) = min(one, fice(i)) Sfcprop%tisfc(i) = tice(i) else ! this would be over open ocean or land (no ice fraction) Sfcprop%hice(i) = zero From fdb78aedadcb44b5263f753c70817fad10e3fc58 Mon Sep 17 00:00:00 2001 From: "Shan.Sun" Date: Sat, 11 Jan 2020 00:23:03 +0000 Subject: [PATCH 04/15] Cosmetic changes --- atmos_model.F90 | 2 +- gfsphysics/GFS_layer/GFS_physics_driver.F90 | 4 ++-- io/FV3GFS_io.F90 | 12 +++++------- 3 files changed, 8 insertions(+), 10 deletions(-) diff --git a/atmos_model.F90 b/atmos_model.F90 index af24b1511..55d0c7974 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -1697,7 +1697,7 @@ subroutine assign_importdata(rc) ix = Atm_block%ixp(i,j) IPD_Data(nb)%Coupling%ficein_cpl(ix) = zero if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then - datar8(i,j) = datar8(i,j)/IPD_Data(nb)%Sfcprop%oceanfrac(ix) !LHS: ice frac wrt non-land portion (not whole cell) + datar8(i,j) = datar8(i,j)/IPD_Data(nb)%Sfcprop%oceanfrac(ix) !convert ice frac wrt whole cell to water area if (datar8(i,j) >= IPD_control%min_seaice) then IPD_Data(nb)%Coupling%ficein_cpl(ix) = datar8(i,j) if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) == one) IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index f40f3162d..d48b5038a 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -1054,9 +1054,9 @@ subroutine GFS_physics_driver & if (slopetyp(i) < 1) slopetyp(i) = 1 endif ! --- ... xw: transfer ice thickness & concentration from global to local variables - zice(i) = Sfcprop%hice(i) Sfcprop%fice(i) = min(one,Sfcprop%fice(i)) ! sea/lake ice fraction wrt water portion, not whole cell fice(i) = Sfcprop%fice(i) + zice(i) = Sfcprop%hice(i) tice(i) = Sfcprop%tisfc(i) ! !GFDL work1(i) = (log(coslat(i) / (nlons(i)*latr)) - dxmin) * dxinv @@ -1986,7 +1986,7 @@ subroutine GFS_physics_driver & if (.not. flag_cice(i)) then if (islmsk(i) == 2) then ! return updated lake ice thickness & concentration to global array Sfcprop%hice(i) = zice(i) - Sfcprop%fice(i) = min(one, fice(i)) + Sfcprop%fice(i) = fice(i) Sfcprop%tisfc(i) = tice(i) else ! this would be over open ocean or land (no ice fraction) Sfcprop%hice(i) = zero diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index 87d23e624..69a7ebbb9 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -1100,23 +1100,21 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) enddo endif - if(Model%cplflx) then + if(Model%cplflx .or. Model%frac_grid) then if (nint(sfc_var2(1,1,33)) == -9999) then - if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing tsfcl') - !--- compute tsfcl from existing variables + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing tsfcl') do nb = 1, Atm_block%nblks do ix = 1, Atm_block%blksz(nb) - Sfcprop(nb)%tsfcl(ix) = Sfcprop(nb)%tsfco(ix) + Sfcprop(nb)%tsfcl(ix) = Sfcprop(nb)%tsfco(ix) !--- compute tsfcl from existing variables enddo enddo endif if (nint(sfc_var2(1,1,34)) == -9999) then - if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorll') - !--- compute zorll from existing variables + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorll') do nb = 1, Atm_block%nblks do ix = 1, Atm_block%blksz(nb) - Sfcprop(nb)%zorll(ix) = Sfcprop(nb)%zorlo(ix) + Sfcprop(nb)%zorll(ix) = Sfcprop(nb)%zorlo(ix) !--- compute zorll from existing variables enddo enddo endif From 8f057b849b9a7495e207a1cf25c44ebb8e27e925 Mon Sep 17 00:00:00 2001 From: "Shan.Sun" Date: Wed, 5 Feb 2020 04:54:24 +0000 Subject: [PATCH 05/15] -- Merged atmos_model.F90, GFS_physics_driver.F90 & FV3GFS_io.F90 to master; -- No changes in results in the stand-alone control test. --- atmos_model.F90 | 40 +- gfsphysics/GFS_layer/GFS_physics_driver.F90 | 506 +++++++++++++++----- io/FV3GFS_io.F90 | 81 +++- 3 files changed, 484 insertions(+), 143 deletions(-) diff --git a/atmos_model.F90 b/atmos_model.F90 index 55d0c7974..e5d39c336 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -113,7 +113,8 @@ module atmos_model_mod FV3GFS_diag_register, FV3GFS_diag_output, & DIAG_SIZE use fv_iau_mod, only: iau_external_data_type,getiauforcing,iau_initialize -use module_fv3_config, only: output_1st_tstep_rst, first_kdt, nsout +use module_fv3_config, only: output_1st_tstep_rst, first_kdt, nsout, & + frestart, restart_endfcst !----------------------------------------------------------------------- @@ -221,7 +222,8 @@ module atmos_model_mod logical,parameter :: flip_vc = .true. #endif - real(kind=IPD_kind_phys), parameter :: zero=0.0, one=1.0 + real(kind=IPD_kind_phys), parameter :: zero = 0.0_IPD_kind_phys, & + one = 1.0_IPD_kind_phys contains @@ -944,7 +946,7 @@ end subroutine update_atmos_model_state subroutine atmos_model_end (Atmos) type (atmos_data_type), intent(inout) :: Atmos !---local variables - integer :: idx + integer :: idx, seconds #ifdef CCPP integer :: ierr #endif @@ -952,9 +954,11 @@ subroutine atmos_model_end (Atmos) !----------------------------------------------------------------------- !---- termination routine for atmospheric model ---- - call atmosphere_end (Atmos % Time, Atmos%grid) - call FV3GFS_restart_write (IPD_Data, IPD_Restart, Atm_block, & - IPD_Control, Atmos%domain) + call atmosphere_end (Atmos % Time, Atmos%grid, restart_endfcst) + if(restart_endfcst) then + call FV3GFS_restart_write (IPD_Data, IPD_Restart, Atm_block, & + IPD_Control, Atmos%domain) + endif #ifdef CCPP ! Fast physics (from dynamics) are finalized in atmosphere_end above; @@ -1457,6 +1461,24 @@ subroutine update_atmos_chemistry(state, rc) enddo enddo + ! -- zero out accumulated fields +!$OMP parallel do default (none) & +!$OMP shared (nj, ni, Atm_block, IPD_Control, IPD_Data) & +!$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)%coupling%rainc_cpl(ix) = zero + if (.not.IPD_Control%cplflx) then + IPD_Data(nb)%coupling%rain_cpl(ix) = zero + IPD_Data(nb)%coupling%snow_cpl(ix) = zero + end if + enddo + enddo + if (IPD_Control%debug) then ! -- diagnostics write(6,'("update_atmos: prsi - min/max/avg",3g16.6)') minval(prsi), maxval(prsi), sum(prsi)/size(prsi) @@ -1697,14 +1719,14 @@ subroutine assign_importdata(rc) ix = Atm_block%ixp(i,j) IPD_Data(nb)%Coupling%ficein_cpl(ix) = zero if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then - datar8(i,j) = datar8(i,j)/IPD_Data(nb)%Sfcprop%oceanfrac(ix) !convert ice frac wrt whole cell to water area - if (datar8(i,j) >= IPD_control%min_seaice) then - IPD_Data(nb)%Coupling%ficein_cpl(ix) = datar8(i,j) + IPD_Data(nb)%Coupling%ficein_cpl(ix) = max(zero, min(one, datar8(i,j)/IPD_Data(nb)%Sfcprop%oceanfrac(ix))) !LHS: ice frac wrt water area + if (IPD_Data(nb)%Coupling%ficein_cpl(ix) >= IPD_control%min_seaice) then if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) == one) IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points IPD_Data(nb)%Coupling%slimskin_cpl(ix) = 4. else if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) == one) IPD_Data(nb)%Sfcprop%slmsk(ix) = zero IPD_Data(nb)%Coupling%slimskin_cpl(ix) = zero + IPD_Data(nb)%Coupling%ficein_cpl(ix) = zero endif endif enddo diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index d48b5038a..e6e0f4c15 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -474,6 +474,13 @@ subroutine GFS_physics_driver & type(GFS_diag_type), intent(inout) :: Diag #endif ! +!## CCPP ## Note: Variables defined locally in this file for temporary calculations +! or transfer of data between schemes are defined in gfsphysics/GFS_layer/GFS_typedefs.F90 +! in the GFS_interstitial_type datatype. Type-bound procedures create, rad_reset, +! phys_reset, and mprint exist to allocate memory, to reset variables used in GFS_radiation_driver.F90, +! to reset variables used in GFS_physics_driver.F90, and to print the contents of the +! data type to the console + ! --- local variables !--- INTEGER VARIABLES @@ -569,10 +576,6 @@ subroutine GFS_physics_driver & real(kind=kind_phys), allocatable, dimension(:,:) :: & savet, saveq, saveu, savev -!--- pass precip type from MP to Noah MP - real(kind=kind_phys), dimension(size(Grid%xlon,1)) :: & - rainn_mp, rainc_mp, snow_mp, graupel_mp, ice_mp - !--- GFDL modification for FV3 real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs+1) ::& @@ -692,22 +695,27 @@ subroutine GFS_physics_driver & !Stateout%gv0(:,:) = Statein%vgrs(:,:) !Stateout%gq0(:,:,:) = Statein%qgrs(:,:,:) +!## CCPP ## Note: Setting local variables from the Model DDT (without additional +! logic attached) is not necessary with the CCPP interstitial schemes with exceptions +! noted below. + !===> ... begin here ldiag_ugwp = Model%ldiag_ugwp ! !===> master = Model%master - me = Model%me - ix = size(Grid%xlon,1) - im = size(Grid%xlon,1) - ipr = min(im,10) + me = Model%me + ix = size(Grid%xlon,1) !## CCPP ## set in GFS_typedefs.F90/interstitial_create + im = size(Grid%xlon,1) !## CCPP ## set in GFS_typedefs.F90/interstitial_create + ipr = min(im,10) !## CCPP ## set in GFS_typedefs.F90/interstitial_create levs = Model%levs lsoil = Model%lsoil ntrac = Model%ntrac dtf = Model%dtf dtp = Model%dtp +!## CCPP ##* this block not yet in CCPP !------- ! For COORDE-2019 averaging with fwindow, it was done before ! 3Diag fixes and averaging ingested using "fdaily"-factor @@ -722,9 +730,11 @@ subroutine GFS_physics_driver & print *, 'VAY Model%fhzero = 0., Bad Averaged-diagnostics ' endif !------- +!*## CCPP ## kdt = Model%kdt lprnt = Model%lprnt +!## CCPP ## see GFS_typedefs.F90/interstitial_setup_tracers for logic for setting nvdiff nvdiff = ntrac ! vertical diffusion of all tracers! ntcw = Model%ntcw ntoz = Model%ntoz @@ -748,6 +758,7 @@ subroutine GFS_physics_driver & imp_physics = Model%imp_physics +!## CCPP ##* GFS_typedefs.F90/interstitial_setup_tracers nncl = ncld ! perform aerosol convective transport and PBL diffusion @@ -806,11 +817,12 @@ subroutine GFS_physics_driver & if (trans_aero) nvdiff = nvdiff + Model%ntchm if (ntke > 0) nvdiff = nvdiff + 1 ! adding tke to the list endif +!*## CCPP ## ! - ! For CCPP, this is in GFS_Interstitial%phys_reset(Model) in GFS_typedefs.F90 - +!## CCPP ##* GFS_typedefs.F90/interstitial_phys_reset kdtminus1 = kdt - 1 reset = mod(kdtminus1, nint(Model%avg_max_length/dtp)) == 0 +!*## CCPP ## ! !------------------------------------------------------------------------------------------- @@ -847,8 +859,10 @@ subroutine GFS_physics_driver & frain = dtf / dtp +!## CCPP ##* GFS_typedefs.F90/interstitial_create skip_macro = .false. - +!*## CCPP ## +!## CCPP ##* GFS_typedefs.F90/interstitial_setup_tracers if (ntiw > 0) then if (ntclamt > 0) then nn = ntrac - 2 @@ -860,8 +874,11 @@ subroutine GFS_physics_driver & else nn = ntrac + 1 endif +!*## CCPP ## +!## CCPP ##* GFS_typedefs.F90/interstitial_create allocate (clw(ix,levs,nn)) - +!*## CCPP ## +!## CCPP ##* GFS_typedefs.F90/interstitial_create Note: cnvc and cnvw are always allocated and initialized regardless of test condition if (Model%imfdeepcnv >= 0 .or. Model%imfshalcnv > 0 .or. & (Model%npdf3d == 3 .and. Model%num_p3d == 4) .or. & (Model%npdf3d == 0 .and. Model%ncnvcld3d == 1) ) then @@ -872,6 +889,8 @@ subroutine GFS_physics_driver & cnvw(i,k) = zero enddo enddo +!*## CCPP ## +!## CCPP ##* GFS_typedefs.F90/control_initialize Note: these are calculated regardless of test condition if (Model%npdf3d == 3 .and. Model%num_p3d == 4) then num2 = Model%num_p3d + 2 num3 = num2 + 1 @@ -881,32 +900,9 @@ subroutine GFS_physics_driver & !CCPP: num2 = Model%ncnvw !CCPP: num3 = Model%ncnvc endif -! -! --- initization for those precip type used in Noah MP -! - if (Model%lsm == Model%lsm_noahmp) then - do i=1,im - rainn_mp(i) = zero - rainc_mp(i) = zero - snow_mp(i) = zero - graupel_mp(i) = zero - ice_mp(i) = zero - enddo -! --- get the amount of different precip type for Noah MP -! --- convert from m/dtp to mm/s - if (Model%imp_physics == Model%imp_physics_mg .or. & - Model%imp_physics == Model%imp_physics_gfdl) then - tem = one / (dtp*con_p001) - do i=1,im - rainn_mp(i) = tem * (Diag%rain(i)-Diag%rainc(i)) - rainc_mp(i) = tem * Diag%rainc(i) - snow_mp(i) = tem * Diag%snow(i) - graupel_mp(i) = tem * Diag%graupel(i) - ice_mp(i) = tem * Diag%ice(i) - enddo - endif - endif ! if (Model%lsm == Model%lsm_noahmp) +!*## CCPP ## +!## CCPP ##* GFS_surface_generic.F90/GFS_surface_generic_pre_run ! --- set initial quantities for stochastic physics deltas if (Model%do_sppt) then Tbd%dtdtr = zero @@ -919,6 +915,8 @@ subroutine GFS_physics_driver & ! mg, sfc-perts ! --- scale random patterns for surface perturbations with perturbation size ! --- turn vegetation fraction pattern into percentile pattern +!## CCPP ##* Note: initialzations to zero are not needed in GFS_surface_generic.F90/GFS_surface_generic_pre_run +! since this function occurs in GFS_typedefs.F90/interstitial_phys_reset do i=1,im z01d(i) = zero zt1d(i) = zero @@ -956,7 +954,9 @@ subroutine GFS_physics_driver & enddo endif endif +!*## CCPP ## ! +!## CCPP ##* GFS_typedefs.F90/interstitial_create if (Model%do_shoc) then allocate (qrn(im,levs), qsnw(im,levs), & ncpl(im,levs), ncpi(im,levs)) @@ -969,7 +969,7 @@ subroutine GFS_physics_driver & enddo enddo endif -! +!## CCPP ##* GFS_typedefs.F90/coupling_create ## if (imp_physics == Model%imp_physics_thompson) then if(Model%ltaerosol) then allocate(ice00(im,levs)) @@ -979,7 +979,8 @@ subroutine GFS_physics_driver & allocate(ice00(im,levs)) endif endif - +!*## CCPP ## +!## CCPP ##* allocated in GFS_typedefs.F90/interstitial_create; initialized in GFS_typedefs.F90/interstitial_phys_reset if (imp_physics == Model%imp_physics_mg) then ! For MGB double moment microphysics allocate (qlcn(im,levs), qicn(im,levs), w_upi(im,levs), & cf_upi(im,levs), CNV_MFD(im,levs), & @@ -1001,12 +1002,14 @@ subroutine GFS_physics_driver & ncgl(i,k) = zero enddo enddo -! +!*## CCPP ## +!## CCPP ##* These variables are currently being allocated fully (im,levs) in GFS_typedefs.F90/interstitial_create else allocate (qlcn(1,1), qicn(1,1), w_upi(1,1), cf_upi(1,1), & CNV_MFD(1,1), CNV_DQLDT(1,1), & ! CNV_MFD(1,1), CNV_PRC3(1,1), CNV_DQLDT(1,1), & clcn(1,1), cnv_fice(1,1), cnv_ndrop(1,1), cnv_nice(1,1)) +!## CCPP ##* The following variables are local to gfdl_cloud_microphys.F90/gfdl_cloud_microphys_run if (imp_physics == Model%imp_physics_gfdl) then ! GFDL MP allocate (delp(im,1,levs), dz(im,1,levs), uin(im,1,levs), & vin(im,1,levs), pt(im,1,levs), qv1(im,1,levs), ql1(im,1,levs), & @@ -1017,7 +1020,9 @@ subroutine GFS_physics_driver & qg_dt(im,1,levs), p123(im,1,levs), refl(im,1,levs), den(im,levs)) endif endif +!*## CCPP ## +!## CCPP ## Only get_prs_fv3.F90/get_prs_fv3_run is a scheme (GFS_HYDRO is assumed to be undefined) #ifdef GFS_HYDRO call get_prs(im, ix, levs, ntrac, Statein%tgrs, Statein%qgrs, & Model%thermodyn_id, Model%sfcpress_id, & @@ -1028,7 +1033,9 @@ subroutine GFS_physics_driver & call get_prs_fv3 (ix, levs, ntrac, Statein%phii, Statein%prsi, & Statein%tgrs, Statein%qgrs, del, del_gz) #endif +!*## CCPP ## +!## CCPP ##* GFS_surface_generic.F90/GFS_surface_generic_pre_run do i = 1, IM sigmaf(i) = max( Sfcprop%vfrac(i),0.01 ) islmsk(i) = nint(Sfcprop%slmsk(i)) @@ -1053,10 +1060,14 @@ subroutine GFS_physics_driver & if (vegtype(i) < 1) vegtype(i) = 17 if (slopetyp(i) < 1) slopetyp(i) = 1 endif +!*## CCPP ## ! --- ... xw: transfer ice thickness & concentration from global to local variables +!## CCPP ## global to local variable transfer not necessary for these two Sfcprop%fice(i) = min(one,Sfcprop%fice(i)) ! sea/lake ice fraction wrt water portion, not whole cell fice(i) = Sfcprop%fice(i) zice(i) = Sfcprop%hice(i) +!*## CCPP ##* +!## CCPP ##* GFS_surface_composites.F90/GFS_surface_composites_pre_run tice(i) = Sfcprop%tisfc(i) ! !GFDL work1(i) = (log(coslat(i) / (nlons(i)*latr)) - dxmin) * dxinv @@ -1066,23 +1077,34 @@ subroutine GFS_physics_driver & work1(i) = max(zero, min(one, work1(i))) work2(i) = one - work1(i) Diag%psurf(i) = Statein%pgr(i) +!*## CCPP ## +!## CCPP ##* GFS_surface_generic.F90/GFS_surface_generic_pre_run work3(i) = Statein%prsik(i,1) / Statein%prslk(i,1) +!*## CCPP ## !GFDL tem1 = con_rerth * (con_pi+con_pi)*coslat(i)/nlons(i) !GFDL tem2 = con_rerth * con_pi / latr !GFDL garea(i) = tem1 * tem2 +!## CCPP ## global to local variable transfer not necessary for these variables tem1 = Grid%dx(i) tem2 = Grid%dx(i) garea(i) = Grid%area(i) +!*## CCPP ## +!## CCPP ##* gwdc.f/gwdc_pre_run dlength(i) = sqrt( tem1*tem1+tem2*tem2 ) cldf(i) = Model%cgwf(1) * work1(i) + Model%cgwf(2) * work2(i) +!*## CCPP ## +!## CCPP ##* cs_conv.F90/cs_conv_pre_run wcbmax(i) = Model%cs_parm(1) * work1(i) + Model%cs_parm(2) * work2(i) -! +!*## CCPP ## +!## CCPP ##* GFS_typedefs.F90/interstitial_phys_reset dry(i) = .false. icy(i) = .false. wet(i) = .false. flag_cice(i) = .false. +!*## CCPP ## enddo ! +!## CCPP ##* GFS_surface_generic.F90/GFS_surface_generic_pre_run if (Model%cplflx) then do i=1,im islmsk_cice(i) = nint(Coupling%slimskin_cpl(i)) @@ -1097,8 +1119,9 @@ subroutine GFS_physics_driver & endif enddo endif +!*## CCPP ## -! DH* In CCPP, this is in GFS_surface_composites_pre +!## CCPP ##* GFS_surface_composites.F90/GFS_surface_composites_pre if (Model%frac_grid) then do i = 1, IM frland(i) = Sfcprop%landfrac(i) @@ -1117,11 +1140,13 @@ subroutine GFS_physics_driver & fice(i) = zero endif endif - if (icy(i)) Sfcprop%tsfco(i) = max(Sfcprop%tsfco(i), Sfcprop%tisfc(i), tgice) + if (fice(i) < one ) then + wet(i)=.true. !there is some open ocean/lake water! + if (.not. Model%cplflx) Sfcprop%tsfco(i) = max(Sfcprop%tsfco(i), Sfcprop%tisfc(i), tgice) + end if else fice(i) = zero endif - if (frland(i)+fice(i)*(one-frland(i)) < one) wet(i)=.true. !there is some open ocean/lake water! enddo else do i = 1, IM @@ -1181,7 +1206,6 @@ subroutine GFS_physics_driver & enddo enddo -! DH* In CCPP, this is in GFS_surface_composites_pre if (.not. Model%cplflx .or. .not. Model%frac_grid) then do i=1,im Sfcprop%zorll(i) = Sfcprop%zorl(i) @@ -1223,9 +1247,9 @@ subroutine GFS_physics_driver & semis3(i,2) = 0.95d0 endif enddo -! *DH +!*## CCPP ## -! +!## CCPP ## global to local variable transfer not necessary for these variables ! --- ... transfer soil moisture and temperature from global to local variables do k=1,lsoil do i=1,im @@ -1234,6 +1258,7 @@ subroutine GFS_physics_driver & slsoil(i,k) = Sfcprop%slc(i,k) !! clu: slc -> slsoil enddo enddo +!*## CCPP ## do k=1,levs do i=1,im @@ -1242,6 +1267,7 @@ subroutine GFS_physics_driver & dtdt(i,k) = zero dtdtc(i,k) = zero +!## CCPP ##* GFS_typedefs.F90/interstitial_phys_reset !vay-2018 ! Pure tendency arrays w/o accumulation of Phys-tendencies from each ! chain of GFS-physics (later add container for species) @@ -1257,9 +1283,10 @@ subroutine GFS_physics_driver & gw_dvdt(i,k) = zero gw_dtdt(i,k) = zero gw_kdis(i,k) = zero -! *DH +!*## CCPP ## enddo enddo +!## CCPP ##* GFS_suite_interstitial.F90/GFS_suite_interstitial_1_run do n=1,ntrac do k=1,levs do i=1,im @@ -1267,7 +1294,9 @@ subroutine GFS_physics_driver & enddo enddo enddo +!*## CCPP ## +!## CCPP ##* This block is not yet in CCPP. !----------------------------------------------- !vay-2018-19 ORO/UGWP process-oriented diagnostics ! @@ -1373,6 +1402,7 @@ subroutine GFS_physics_driver & endif endif !===========================Above Phys-tend Diag for COORDE ====================== +!*## CCPP ## ! --- ... initialize dtdt with heating rate from dcyc2 @@ -1380,7 +1410,7 @@ subroutine GFS_physics_driver & ! faster model time steps. ! sw: using cos of zenith angle as scaling factor ! lw: using surface air skin temperature as scaling factor - +!## CCPP ##* This is not in the CCPP yet. if (Model%pre_rad) then call dcyc2t3_pre_rad & ! --- inputs: @@ -1400,7 +1430,8 @@ subroutine GFS_physics_driver & ) else - +!*## CCPP ## +!** CCPP ## dcyc2.f/dcyc2t3_run Note: Check for Model%pre_rad was omitted, so this option is broken in CCPP call dcyc2t3 & ! --- inputs: ( Model%solhr, Model%slag, Model%sdec, Model%cdec, Grid%sinlat, & @@ -1421,15 +1452,17 @@ subroutine GFS_physics_driver & adjnirbmu, adjnirdfu, adjvisbmu, adjvisdfu, & adjnirbmd, adjnirdfd, adjvisbmd, adjvisdfd & ) - +!*## CCPP ## ! ! save temp change due to radiation - need for sttp stochastic physics !--------------------------------------------------------------------- endif ! +!## CCPP ##* This is not in the CCPP yet. if (Model%lsidea) then !idea jw dtdt(:,:) = zero endif +!*## CCPP ## ! --- convert lw fluxes for land/ocean/sea-ice models ! note: for sw: adjsfcdsw and adjsfcnsw are zenith angle adjusted downward/net fluxes. @@ -1449,12 +1482,15 @@ subroutine GFS_physics_driver & ! --- ... define the downward lw flux absorbed by ground +!## CCPP ##* GFS_surface_composites.F90/GFS_surface_composites_pre_run do i=1,im if (dry(i)) gabsbdlw3(i,1) = semis3(i,1) * adjsfcdlw(i) if (icy(i)) gabsbdlw3(i,2) = semis3(i,2) * adjsfcdlw(i) if (wet(i)) gabsbdlw3(i,3) = semis3(i,3) * adjsfcdlw(i) enddo +!*## CCPP ## +!## CCPP ##* GFS_suite_interstitial.F90/GFS_suite_interstitial_2_run if (Model%lssav) then ! --- ... accumulate/save output variables ! --- ... sunshine duration time is defined as the length of time (in mdl output @@ -1542,8 +1578,8 @@ subroutine GFS_physics_driver & endif ! end if_lssav_block do i=1,im - kcnv(i) = 0 - kinver(i) = levs + kcnv(i) = 0 !## CCPP ## GFS_typedefs.F90/interstitial_phys_reset + kinver(i) = levs !## CCPP ## GFS_typedefs.F90/interstitial_phys_reset invrsn(i) = .false. tx1(i) = zero tx2(i) = 10.0 @@ -1592,10 +1628,11 @@ subroutine GFS_physics_driver & enddo enddo endif - +!*## CCPP ## ! --- ... lu: initialize flag_guess, flag_iter, tsurf +!## CCPP ##* These initializations are done in GFS_typedefs.F90/interstitial_phys_reset except for as noted below do i=1,im ! tsurf(i) = Sfcprop%tsfc(i) flag_guess(i) = .false. @@ -1620,16 +1657,19 @@ subroutine GFS_physics_driver & Statein%vgrs(i,1)*Statein%vgrs(i,1)) & + max(zero, min(Tbd%phy_f2d(i,Model%num_p2d), 30.0)), one) enddo +!*## CCPP ## ! --- ... lu: iter-loop over (sfc_diff,sfc_drv,sfc_ocean,sfc_sice) - +!## CCPP ##* This loop is implemented using the subcycle/iteration capability in the CCPP SDF do iter=1,2 +!*## CCPP ## ! --- ... surface exchange coefficients ! ! if (lprnt) write(0,*)' tsfc=',Sfcprop%tsfc(ipr),' tsurf=',tsurf(ipr),'iter=', & ! iter ,'wet=',wet(ipr),'dry=',dry(ipr),' icy=',icy(ipr) +!## CCPP ##* sfc_diff.f/sfc_diff_run call sfc_diff & ! --- inputs: (im, Statein%pgr, & @@ -1647,13 +1687,16 @@ subroutine GFS_physics_driver & ! cd3, cdq3, rb3, stress3, ffmm3, ffhh3, fm103, fh23, wind, lprnt, ipr) ! ! --- ... lu: update flag_guess - +!*## CCPP ## +!## CCPP ##* GFS_surface_loop_control/GFS_surface_loop_control_part1_run do i=1,im if (iter == 1 .and. wind(i) < 2.0) then flag_guess(i) = .true. endif enddo - +!*## CCPP ## +!## CCPP ##* sfc_nst.f/sfc_nst_pre_run Note: the conditional is not included in the CCPP scheme, so calling +! this code is controlled by its presence in the active CCPP SDF if (Model%nstf_name(1) > 0) then do i=1,im if (wet(i)) then @@ -1684,7 +1727,8 @@ subroutine GFS_physics_driver & endif ! if (lprnt) write(0,*)' bef nst tseal=',tseal(ipr) & ! ,' tsfc3=',tsfc3(ipr,3),' tsurf3=',tsurf3(ipr,3),' tem=',tem - +!*## CCPP ## +!## CCPP ##* sfc_nst.f/sfc_nst_run call sfc_nst & ! --- inputs: (im, Statein%pgr, Statein%ugrs(:,1), Statein%vgrs(:,1), & @@ -1705,6 +1749,8 @@ subroutine GFS_physics_driver & ! --- outputs: qss3(:,3), gflx3(:,3), cmm3(:,3), chh3(:,3), evap3(:,3), & hflx3(:,3), ep1d3(:,3)) +!*## CCPP ## +!## CCPP ##* sfc_nst.f/sfc_nst_post_run ! do i=1,im !! if (wet(i) .and. .not.icy(i)) then @@ -1733,14 +1779,15 @@ subroutine GFS_physics_driver & endif enddo endif - +!*## CCPP ## ! if (lprnt) print *,' tseaz2=',Sfcprop%tsfc(ipr),' tref=',tref(ipr), & ! & ' dt_cool=',dt_cool(ipr),' dt_warm=',dt_warm(ipr),' kdt=',kdt +!## CCPP ## Note: This conditional is replaced by whether the sfc_ocean scheme is in the CCPP SDF else ! --- ... surface energy balance over ocean - +!## CCPP ##* sfc_ocean.F/sfc_ocean_run call sfc_ocean & ! --- inputs: (im, Statein%pgr, & @@ -1750,6 +1797,7 @@ subroutine GFS_physics_driver & ! --- outputs: qss3(:,3), cmm3(:,3), chh3(:,3), gflx3(:,3), evap3(:,3), & hflx3(:,3), ep1d3(:,3)) +!*## CCPP ## endif ! if nstf_name(1) > 0 @@ -1760,12 +1808,16 @@ subroutine GFS_physics_driver & ! --- ... surface energy balance over land ! +!## CCPP ##* Note: the conditional is not included in the CCPP, so calling +! the LSM scheme is controlled by its presence in the active CCPP SDF if (Model%lsm == Model%lsm_noah) then ! noah lsm call +!*## CCPP ## ! if (lprnt) write(0,*)' tseal=',tseal(ipr),' tsurf=',tsurf(ipr),iter & ! ,' stsoil0=',stsoil(ipr,:) ! &,' pgr=',pgr(ipr),' sfcemis=',sfcemis(ipr) +!## CCPP ##* sfc_drv.f/lsm_noah_run call sfc_drv & ! --- inputs: (im, lsoil, Statein%pgr, & @@ -1787,12 +1839,14 @@ subroutine GFS_physics_driver & hflx3(:,1), ep1d3(:,1), runof, & cmm3(:,1), chh3(:,1), evbs, evcw, sbsno, snowc, Diag%soilm, & snohf, Diag%smcwlt2, Diag%smcref2, Diag%wet1) +!*## CCPP ## ! if (lprnt) write(0,*)' tseae=',tseal(ipr),' tsurf=',tsurf(ipr),iter& ! ,' phy_f2d=',phy_f2d(ipr,num_p2d) ! if (lprnt) write(0,*)' hflx3=',hflx3(ipr,:),' evap3=',evap3(i,:) +!## CCPP ##* sfc_noahmp_drv.f/noahmpdrv_run ! Noah MP call ! elseif (Model%lsm == Model%lsm_noahmp) then @@ -1812,8 +1866,8 @@ subroutine GFS_physics_driver & Model%iopt_inf, Model%iopt_rad, Model%iopt_alb, & Model%iopt_snf, Model%iopt_tbot, Model%iopt_stc, & grid%xlat, xcosz, Model%yearlen, Model%julian, Model%imn,& - rainn_mp, rainc_mp, snow_mp, graupel_mp, ice_mp, & - + Sfcprop%drainncprv, Sfcprop%draincprv, Sfcprop%dsnowprv, & + Sfcprop%dgraupelprv, Sfcprop%diceprv, & ! --- in/outs: weasd3(:,1), snowd3(:,1), tsfc3(:,1), tprcp3(:,1), & Sfcprop%srflag, smsoil, stsoil, slsoil, Sfcprop%canopy, & @@ -1838,6 +1892,7 @@ subroutine GFS_physics_driver & ! if (lprnt) write(0,*)' tseae=',tsea(ipr),' tsurf=',tsurf(ipr),iter & ! &,' phy_f2d=',phy_f2d(ipr,num_p2d) +!*## CCPP ## elseif (Model%lsm == Model%lsm_ruc) then write (0,*) 'RUC LSM is available only in CCPP' @@ -1850,14 +1905,16 @@ subroutine GFS_physics_driver & ! &,' stsoil=',stsoil(ipr,:) ! --- ... surface energy balance over seaice - +!## CCPP ##* sfc_sice.f/sfc_sice_run (local adjustment to avoid resetting islmsk after call to sfc_sice_run) if (Model%cplflx) then do i=1,im if (flag_cice(i)) then islmsk (i) = islmsk_cice(i) endif enddo +!*## CCPP ## +!## CCPP ##* sfc_cice.f/sfc_cice_run ! call sfc_cice for sea ice points in the coupled model (i.e. islmsk=4) ! call sfc_cice & @@ -1871,10 +1928,12 @@ subroutine GFS_physics_driver & qss3(:,2), cmm3(:,2), chh3(:,2), evap3(:,2), hflx3(:,2), & stress3(:,2)) endif +!*## CCPP ## ! ! call sfc_sice for lake ice and for the uncoupled case, sea ice (i.e. islmsk=2) ! +!## CCPP ##* sfc_sice.f/sfc_sice_run call sfc_sice & ! --- inputs: (im, lsoil, Statein%pgr, & @@ -1890,7 +1949,8 @@ subroutine GFS_physics_driver & ! --- outputs: snowd3(:,2), qss3(:,2), snowmt, gflx3(:,2), cmm3(:,2), chh3(:,2), & evap3(:,2), hflx3(:,2)) - +!*## CCPP ## +!## CCPP ##* This section is not needed for CCPP. if (Model%cplflx) then do i = 1, im if (flag_cice(i)) then @@ -1898,12 +1958,13 @@ subroutine GFS_physics_driver & endif enddo endif +!*## CCPP ## ! if (lprnt) write(0,*)' tseaafticemodel =',tsfc3(ipr,2),' me=',me & ! &, ' kdt=',kdt,' iter=',iter,' fice=',fice(ipr) ! --- ... lu: update flag_iter and flag_guess - +!## CCPP ##* GFS_surface_loop_control.F90/GFS_surface_loop_control_part_2 do i=1,im flag_iter(i) = .false. flag_guess(i) = .false. @@ -1916,12 +1977,14 @@ subroutine GFS_physics_driver & endif enddo +!*## CCPP ## enddo ! end iter_loop ! --- generate ocean/land/ice composites +!## CCPP ##* GFS_surface_compoistes.F90/GFS_surface_composites_post_run if (Model%frac_grid) then do i=1, im ! @@ -2044,7 +2107,7 @@ subroutine GFS_physics_driver & Sfcprop%zorlo(i) = zorl3(i,3) if (flag_cice(i)) then ! this was already done for lake ice in sfc_sice - txi = Sfcprop%fice(i) + txi = fice(i) txo = one - txi evap(i) = txi * evap3(i,2) + txo * evap3(i,3) hflx(i) = txi * hflx3(i,2) + txo * hflx3(i,3) @@ -2063,32 +2126,43 @@ subroutine GFS_physics_driver & endif enddo endif ! if (Model%frac_grid) +!*## CCPP ## ! --- compositing done ! if (lprnt) write(0,*) 'tisfc=',Sfcprop%tisfc(ipr),'tice=',tice(ipr),' kdt=',kdt do i=1,im +!## CCPP ##* GFS_surface_generic.F90/GFS_surface_generic_post_run Diag%epi(i) = ep1d(i) +!*## CCPP ## Diag%dlwsfci(i) = adjsfcdlw(i) Diag%ulwsfci(i) = adjsfculw(i) +!## CCPP ##* dcyc2.f/dcyc2t3_post_run Diag%uswsfci(i) = adjsfcdsw(i) - adjsfcnsw(i) +!*## CCPP ## Diag%dswsfci(i) = adjsfcdsw(i) +!## CCPP ##* GFS_surface_generic.F90/GFS_surface_generic_post_run Diag%gfluxi(i) = gflx(i) Diag%t1(i) = Statein%tgrs(i,1) Diag%q1(i) = Statein%qgrs(i,1,1) Diag%u1(i) = Statein%ugrs(i,1) Diag%v1(i) = Statein%vgrs(i,1) +!*## CCPP ## enddo ! --- ... update near surface fields +!## CCPP ##* sfc_diag.f/sfc_diag_run call sfc_diag (im, Statein%pgr, Statein%ugrs(:,1), Statein%vgrs(:,1), & Statein%tgrs(:,1), Statein%qgrs(:,1,1), work3, evap, & Sfcprop%ffmm, Sfcprop%ffhh, fm10, fh2, Sfcprop%tsfc, qss, & Sfcprop%f10m, Diag%u10m, Diag%v10m, Sfcprop%t2m, Sfcprop%q2m) +!*## CCPP ## +!## CCPP ##* This block is not in the CCPP Tbd%phy_f2d(:,Model%num_p2d) = zero +!*## CCPP ## if (Model%lsm == Model%lsm_noahmp) then do i=1,im @@ -2177,7 +2251,8 @@ subroutine GFS_physics_driver & Diag%snowca(i) = Diag%snowca(i) + snowc(i) * dtf Diag%snohfa(i) = Diag%snohfa(i) + snohf(i) * dtf Diag%ep(i) = Diag%ep(i) + ep1d(i) * dtf - +!*## CCPP ## +!## CCPP ##* sfc_diag_post.F90/sfc_diag_post_run Diag%tmpmax(i) = max(Diag%tmpmax(i), Sfcprop%t2m(i)) Diag%tmpmin(i) = min(Diag%tmpmin(i), Sfcprop%t2m(i)) @@ -2200,6 +2275,7 @@ subroutine GFS_physics_driver & enddo endif +!*## CCPP ## !!!!!!!!!!!!!!!!!Commented by Moorthi on July 18, 2012 !!!!!!!!!!!!!!!!!!! ! do i=1,im @@ -2238,6 +2314,8 @@ subroutine GFS_physics_driver & ! write(0,*)' before monsho hflx=',hflx,' me=',me ! write(0,*)' before monsho evap=',evap,' me=',me +!## CCPP ##* Note: In the CCPP, the vdftra array is prepared in GFS_PBL_generic.F90/GFS_PBL_generic_pre_run +! regardless of the following conditions. Therefore, this block is redundant in the CCPP and is not included. if (nvdiff == ntrac .or. Model%do_ysu .or. Model%shinhong) then ! ntiwx = 0 @@ -2262,6 +2340,7 @@ subroutine GFS_physics_driver & else if (Model%satmedmf) then if (Model%isatmedmf == 0) then ! initial version of satmedmfvdif (Nov 2018) +!## CCPP ##* satmedmfvdif.F/satmedmfvdif_run Note: The conditional above is checked in satmedmfvdif_init call satmedmfvdif(ix, im, levs, nvdiff, ntcw, ntiw, ntke, & dvdt, dudt, dtdt, dqdt, & Statein%ugrs, Statein%vgrs, Statein%tgrs, Statein%qgrs, & @@ -2272,7 +2351,9 @@ subroutine GFS_physics_driver & Statein%prslk, Statein%phii, Statein%phil, dtp, & Model%dspheat, dusfc1, dvsfc1, dtsfc1, dqsfc1, Diag%hpbl, & kinver, Model%xkzm_m, Model%xkzm_h, Model%xkzm_s) +!*## CCPP ## elseif (Model%isatmedmf == 1) then ! updated version of satmedmfvdif (May 2019) +!## CCPP ##* satmedmfvdifq.F/satmedmfvdifq_run Note: The conditional above is checked in satmedmfvdifq_init call satmedmfvdifq(ix, im, levs, nvdiff, ntcw, ntiw, ntke, & dvdt, dudt, dtdt, dqdt, & Statein%ugrs, Statein%vgrs, Statein%tgrs, Statein%qgrs, & @@ -2284,6 +2365,7 @@ subroutine GFS_physics_driver & Model%dspheat, dusfc1, dvsfc1, dtsfc1, dqsfc1, Diag%hpbl, & kinver, Model%xkzm_m, Model%xkzm_h, Model%xkzm_s, & Model%dspfac, Model%bl_upfr, Model%bl_dnfr) +!*## CCPP ## endif elseif (Model%hybedmf) then if (Model%moninq_fac > 0) then @@ -2357,10 +2439,15 @@ subroutine GFS_physics_driver & endif ! end if_hybedmf endif ! end if_do_shoc else +!*## CCPP ## +!## CCPP ## These variables are allocated in GFS_typedefs.F90/interstitial_create and +! initialized in GFS_typedefs.F90/interstitial_phys_reset; ntiwx is set in +! GFS_typedef.F90/interstitial_setup_tracers allocate(vdftra(ix,levs,nvdiff), dvdftra(im,levs,nvdiff)) dvdftra(:,:,:) = zero ntiwx = 0 ! +!## CCPP ##* GFS_PBL_generic.F90/GFS_PBL_generic_pre_run (ntiwx is set in GFS_typedef.F90/interstitial_setup_tracers) if (imp_physics == Model%imp_physics_wsm6) then ! WSM6 do k=1,levs @@ -2486,10 +2573,12 @@ subroutine GFS_physics_driver & enddo enddo endif -! +!*## CCPP ## ! for SHOC nvdiff=ntrac, so the following is not needed unless cplchm is true ! ----------------------------------------------------- if (Model%do_shoc) then +!## CCPP ##* moninshoc.f/moninshoc_run Note: The conditional above is not checked in the CCPP scheme; +! therefore the use of this scheme is controlled via the CCPP SDF call moninshoc(ix, im, levs, nvdiff, ntcw, nncl, dvdt, dudt, dtdt, dvdftra, & Statein%ugrs, Statein%vgrs, Statein%tgrs, vdftra, & Tbd%phy_f3d(1,1,ntot3d-1), prnum, ntkev, & @@ -2500,9 +2589,11 @@ subroutine GFS_physics_driver & dvsfc1, dtsfc1, dqsfc1, dkt, Diag%hpbl, kinver, & Model%xkzm_m, Model%xkzm_h, Model%xkzm_s, Model%xkzminv, & lprnt, ipr, me) +!*## CCPP ## else if (Model%satmedmf) then if (Model%isatmedmf == 0) then ! initial version of satmedmfvdif (Nov 2018) +!## CCPP ##* satmedmfvdif.F/satmedmfvdif_run Note: The conditional above is checked in satmedmfvdif_init call satmedmfvdif(ix, im, levs, nvdiff, ntcw, ntiwx, ntkev, & dvdt, dudt, dtdt, dvdftra, & Statein%ugrs, Statein%vgrs, Statein%tgrs, vdftra, & @@ -2513,7 +2604,9 @@ subroutine GFS_physics_driver & Statein%prslk, Statein%phii, Statein%phil, dtp, & Model%dspheat, dusfc1, dvsfc1, dtsfc1, dqsfc1, Diag%hpbl, & kinver, Model%xkzm_m, Model%xkzm_h, Model%xkzm_s) +!*## CCPP ## elseif (Model%isatmedmf == 1) then ! updated version of satmedmfvdif (May 2019) +!## CCPP ##* satmedmfvdifq.F/satmedmfvdifq_run Note: The conditional above is checked in satmedmfvdifq_init call satmedmfvdifq(ix, im, levs, nvdiff, ntcw, ntiwx, ntkev, & dvdt, dudt, dtdt, dvdftra, & Statein%ugrs, Statein%vgrs, Statein%tgrs, vdftra, & @@ -2525,8 +2618,11 @@ subroutine GFS_physics_driver & Model%dspheat, dusfc1, dvsfc1, dtsfc1, dqsfc1, Diag%hpbl, & kinver, Model%xkzm_m, Model%xkzm_h, Model%xkzm_s, & Model%dspfac, Model%bl_upfr, Model%bl_dnfr) +!*## CCPP ## endif elseif (Model%hybedmf) then +!## CCPP ## moninedmf.f/hedmf_run Note: The conditional above is not checked in the CCPP scheme; +! therefore the use of this scheme is controlled via the CCPP SDF if ( Model%moninq_fac > 0 ) then call moninedmf(ix, im, levs, nvdiff, ntcw, dvdt, dudt, dtdt, dvdftra, & Statein%ugrs, Statein%vgrs, Statein%tgrs, vdftra, & @@ -2539,6 +2635,8 @@ subroutine GFS_physics_driver & gamt, gamq, dkt, kinver, Model%xkzm_m, Model%xkzm_h, & Model%xkzm_s, lprnt, ipr, & Model%xkzminv, Model%moninq_fac) +!*## CCPP ## +!## CCPP ##* The following schemes are not in the CCPP yet. else call moninedmf_hafs(ix, im, levs, nvdiff, ntcw, dvdt, dudt, dtdt, dvdftra, & Statein%ugrs, Statein%vgrs, Statein%tgrs, vdftra, & @@ -2585,7 +2683,8 @@ subroutine GFS_physics_driver & endif ! end if_satmedmf endif ! end if_do_shoc -! +!*## CCPP ## +!## CCPP ## GFS_PBL_generic.F90/GFS_PBL_generic_post_run if (ntke > 0) then do k=1,levs do i=1,im @@ -2699,7 +2798,9 @@ subroutine GFS_physics_driver & deallocate(vdftra, dvdftra) endif +!*## CCPP ## +!## CCPP ##* GFS_PBL_generic.F90/GFS_PBL_generic_post_run if (Model%cplchm) then do i = 1, im tem1 = max(Diag%q1(i), 1.e-8) @@ -2761,7 +2862,9 @@ subroutine GFS_physics_driver & endif ! Ocean only, NO LAKES enddo endif +!*## CCPP ## !-------------------------------------------------------lssav if loop ---------- +!## CCPP ## GFS_PBL_generic.F90/GFS_PBL_generic_post_run if (Model%lssav) then do i=1,im Diag%dusfc (i) = Diag%dusfc(i) + dusfc1(i)*dtf @@ -2800,7 +2903,9 @@ subroutine GFS_physics_driver & endif endif ! end if_lssav +!*## CCPP ## +!## CCPP ##* This block not yet in CCPP. ! if (ldiag_ugwp) then ! @@ -2827,12 +2932,14 @@ subroutine GFS_physics_driver & enddo endif endif +!*## CCPP ## !============================================================= GW-physics start ! ! Orographic gravity wave drag parameterization ! --------------------------------------------- +!## CCPP ##* GFS_GWD_generic.F90/GFS_GWD_generic_pre_run if (nmtvr == 14) then ! current operational - as of 2014 do i=1,im ! vay-2018 @@ -2888,7 +2995,9 @@ subroutine GFS_physics_driver & gamma = zero ; sigma = zero ; elvmax = zero endif ! end if_nmtvr +!*## CCPP ## +!## CCPP ##* cires_ugwp.F90/cires_ugwp_run - only V0 is implemented ! !===== UGWP-start: two versions V0 (knob_ugwp_version=0) and V1(knob_ugwp_version=1) ! @@ -3034,6 +3143,8 @@ subroutine GFS_physics_driver & !! enddo !! endif +!## CCPP ##* rayleigh_damp.f/rayleigh_damp_run Note: Conditional IS checked +! within the scheme (returns from scheme if condition is not met) ! Rayleigh damping near the model top if( .not. Model%lsidea .and. Model%ral_ts > zero) then call rayleigh_damp(im, ix, im, levs, dvdt, dudt, dtdt, & @@ -3050,6 +3161,9 @@ subroutine GFS_physics_driver & ! Standard accum-Update before "moist physics" by "PBL + GWP + RF" as in GFS/GSM ! +!## CCPP ##* GFS_suite_interstitial.F90/GFS_suite_stateout_update Note: Terms +! containing gw_* are related to the CIRES UGWP code and are not currently in +! this scheme. do k=1,levs do i=1,im Stateout%gt0(i,k) = Statein%tgrs(i,k) + dtdt(i,k) * dtp @@ -3058,10 +3172,9 @@ subroutine GFS_physics_driver & enddo enddo Stateout%gq0(1:im,:,:) = Statein%qgrs(1:im,:,:) + dqdt(1:im,:,:) * dtp +!*## CCPP ## -! *DH - -! DH* this block not yet in CCPP +!## CCPP ##* This is not in the CCPP yet. !================================================================================ ! above: updates of the state by UGWP oro-GWS and RF-damp ! Diag%tav_ugwp & Diag%uav_ugwp(i,k)-Updated U-T state before moist/micro ! physics @@ -3076,24 +3189,30 @@ subroutine GFS_physics_driver & enddo enddo endif +!*## CCPP ## !================================================================================ ! It is not clear Do we need it, "ideaca_up", having stability check inside UGWP-module - +!## CCPP ##* This is not in the CCPP yet. if (Model%lsidea) then ! idea convective adjustment call ideaca_up(Statein%prsi,Stateout%gt0,ix,im,levs+1) endif +!*## CCPP ## ! --- ... ozone physics if (ntoz > 0 .and. ntrac >= ntoz) then if (oz_coeff > 4) then +!## CCPP ##* ozphys_2015.f/ozphys_2015_run Note: The conditionals above are not +! checked in the scheme. The scheme's use is controlled by its presense in the +! CCPP SDF call ozphys_2015 (ix, im, levs, levozp, dtp, & Stateout%gq0(1,1,ntoz), & Stateout%gq0(1,1,ntoz), & Stateout%gt0, oz_pres, Statein%prsl, & Tbd%ozpl, oz_coeff, del, Model%ldiag3d, & dq3dt_loc(1,1,6), me) +!*## CCPP ## ! if (Model%ldiag3d) then ! do k=1,levs ! do i=1,im @@ -3105,12 +3224,14 @@ subroutine GFS_physics_driver & ! enddo ! endif else +!## CCPP ##* ozphys.f/ozphys_run call ozphys (ix, im, levs, levozp, dtp, & Stateout%gq0(1,1,ntoz), & Stateout%gq0(1,1,ntoz), & Stateout%gt0, oz_pres, Statein%prsl, & Tbd%ozpl, oz_coeff, del, Model%ldiag3d, & dq3dt_loc(1,1,6), me) +!*## CCPP ## ! if (Model%ldiag3d) then ! do k=1,levs ! do i=1,im @@ -3125,10 +3246,13 @@ subroutine GFS_physics_driver & endif if (Model%h2o_phys) then +!## CCPP ## h2ophys.f/h2ophys_run Note: The conditional is not checked within +! the scheme. The scheme's use is controlled via the CCPP SDF. call h2ophys (ix, im, levs, levh2o, dtp, Stateout%gq0(1,1,1), & Stateout%gq0(1,1,1), h2o_pres, Statein%prsl, & Tbd%h2opl, h2o_coeff, Model%ldiag3d, & dq3dt_loc(1,1,1), me) +!*## CCPP ## endif ! --- ... to side-step the ozone physics @@ -3163,6 +3287,7 @@ subroutine GFS_physics_driver & ! &,' lat=',lat,' kdt=',kdt,' me=',me ! if (lprnt) write(7000,*)' bef convection gv0=',gv0(ipr,:) +!## CCPP ## GFS_DCNV_generic.F90/GFS_DCNV_generic_pre_run if (Model%ldiag3d) then do k=1,levs do i=1,im @@ -3175,10 +3300,16 @@ subroutine GFS_physics_driver & dtdt(1:im,:) = Stateout%gt0(1:im,:) endif ! end if_ldiag3d/cnvgwd - if (Model%ldiag3d) then + if (Model%ldiag3d .or. Model%cplchm) then dqdt(1:im,:,1) = Stateout%gq0(1:im,:,1) - endif ! end if_ldiag3d + endif ! end if_ldiag3d/cplchm + + if (Model%cplchm) then + Coupling%dqdti(1:im,:) = zero + endif ! end if_cplchm +!*## CCPP ## +!## CCPP ## Only get_prs_fv3.F90/get_phi_fv3_run is a scheme (GFS_HYDRO is assumed to be undefined) #ifdef GFS_HYDRO call get_phi(im, ix, levs, ntrac, Stateout%gt0, Stateout%gq0, & Model%thermodyn_id, Model%sfcpress_id, & @@ -3189,7 +3320,10 @@ subroutine GFS_physics_driver & call get_phi_fv3 (ix, levs, ntrac, Stateout%gt0, Stateout%gq0, & del_gz, Statein%phii, Statein%phil) #endif +!*## CCPP ## +!## CCPP ## These variables are initialized every physics time step through +! GFS_typedefs.F90/interstitial_phys_reset do k=1,levs do i=1,im clw(i,k,1) = zero @@ -3205,16 +3339,20 @@ subroutine GFS_physics_driver & ice00 (:,:) = zero endif endif - +!*## CCPP ## ! --- ... for convective tracer transport (while using ras, csaw, or samf) ! (the code here implicitly assumes that ntiw=ntcw+1) +!## CCPP ## Most of this code block is in GFS_typedefs.F90/interstitial_setup_tracers except +! for code that needs to be executed every time step (noted below). For those lines, +! they are in GFS_suite_interstitial.F90/GFS_suite_interstitial_3_run. ntk = 0 tottracer = 0 if (Model%cscnv .or. Model%satmedmf .or. Model%trans_trac ) then otspt(:,:) = .true. ! otspt is used only for cscnv otspt(1:3,:) = .false. ! this is for sp.hum, ice and liquid water +!## CCPP ##* GFS_suite_interstitial.F90/GFS_suite_interstitial_3_run tracers = 2 do n=2,ntrac if ( n /= ntcw .and. n /= ntiw .and. n /= ntclamt .and. & @@ -3226,6 +3364,7 @@ subroutine GFS_physics_driver & clw(i,k,tracers) = Stateout%gq0(i,k,n) enddo enddo +!*## CCPP ## if (ntke == n ) then otspt(tracers+1,1) = .false. ntk = tracers @@ -3239,19 +3378,21 @@ subroutine GFS_physics_driver & enddo tottracer = tracers - 2 endif ! end if_ras or cfscnv or samf +!*## CCPP ## ! if (kdt == 1 .and. me == 0) & ! write(0,*)' trans_trac=',Model%trans_trac,' tottracer=', & ! & tottracer,' kdt=',kdt,' ntk=',ntk - +!## CCPP ##* These variables are initialized in GFS_typedefs.F90/interstitial_phys_reset do i=1,im ktop(i) = 1 kbot(i) = levs enddo +!*## CCPP ## ! --- ... calling condensation/precipitation processes ! -------------------------------------------- - +!## CCPP ## GFS_suite_interstitial.F90/GFS_suite_interstitial_3_run if (ntcw > 0) then ! if (imp_physics == Model%imp_physics_mg .and. .not. Model%do_shoc) then ! compute rhc for GMAO macro physics cloud pdf if (imp_physics == Model%imp_physics_mg .and. Model%crtrh(2) < 0.5) then ! compute rhc for GMAO macro physics cloud pdf @@ -3298,13 +3439,17 @@ subroutine GFS_physics_driver & enddo endif endif ! ntcw > 0 +!*## CCPP ## ! if (imp_physics == Model%imp_physics_zhao_carr .or. & imp_physics == Model%imp_physics_zhao_carr_pdf) then ! zhao-carr microphysics +!## CCPP ##* precpd.f/zhaocarr_precpd_run do i=1,im psautco_l(i) = Model%psautco(1)*work1(i) + Model%psautco(2)*work2(i) prautco_l(i) = Model%prautco(1)*work1(i) + Model%prautco(2)*work2(i) enddo +!*## CCPP ## +!## CCPP ##* GFS_suite_interstitial.F90/GFS_suite_interstitial_3_run do k=1,levs do i=1,im clw(i,k,1) = Stateout%gq0(i,k,ntcw) @@ -3333,24 +3478,39 @@ subroutine GFS_physics_driver & clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water enddo enddo +!*## CCPP ## +!## CCPP ## These lines are not in the CCPP since it appeared that they were +! not needed. These variables are only ever used if (imp_physics == 99 .or. imp_physics == 98) +! which is handled by the first if statement. else do i=1,im psautco_l(i) = Model%psautco(1)*work1(i) + Model%psautco(2)*work2(i) prautco_l(i) = Model%prautco(1)*work1(i) + Model%prautco(2)*work2(i) enddo +!*## CCPP ## +!## CCPP ##* GFS_suite_interstitial.F90/GFS_suite_interstitial_3_run rhc(:,:) = one +!*## CCPP ## endif ! ! Call SHOC if do_shoc is true and shocaftcnv is false ! +!## CCPP ##* gcm_shoc.F90/shoc_run Note: do_shoc is not checked in the scheme, so +! using this scheme is controlled via the CCPP SDF. if (Model%do_shoc .and. .not. Model%shocaftcnv) then if (imp_physics == Model%imp_physics_mg) then do k=1,levs do i=1,im +!## CCPP ##* These lines are commented out in gcm_shoc.F90/shoc_run since they are +! previously executed in GFS_suite_interstitial.F90/GFS_suite_interstitial_3_run ! clw(i,k,1) = Stateout%gq0(i,k,ntiw) ! ice ! clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water +!*## CCPP ## +!## CCPP ##* These lines are commented out in gcm_shoc.F90/shoc_run since it is +! not necessary to copy global variables to local variables ncpl(i,k) = Stateout%gq0(i,k,ntlnc) ncpi(i,k) = Stateout%gq0(i,k,ntinc) +!*## CCPP ## enddo enddo if (abs(Model%fprcp) == 1 .or. mg3_as_mg2) then @@ -3462,6 +3622,7 @@ subroutine GFS_physics_driver & ! &, Stateout%gq0(1:ix,1:levs,1),clw(1,1,2),clw(1,1,1) & ! &, ' shoc ', grid%xlon(1:im), grid%xlat(1:im)) +!## CCPP ## this is in CCPP's gcm_shoc (but commented out because not needed) if (imp_physics == Model%imp_physics_mg) then do k=1,levs do i=1,im @@ -3470,7 +3631,7 @@ subroutine GFS_physics_driver & enddo enddo endif - +!*## CCPP ## ! do k=1,levs ! do i=1,im ! sgs_cld(i,k) = sgs_cld(i,k) + shoc_cld(i,k) @@ -3493,6 +3654,7 @@ subroutine GFS_physics_driver & ! write(0,*)' aft shoc gq0=',gq0(1,:,1),' lat=',lat ! write(0,*)' aft shoc gu0=',gu0(1,:),' lat=',lat ! +!*## CCPP ## endif ! if(do_shoc) ! @@ -3500,6 +3662,9 @@ subroutine GFS_physics_driver & ! ----------------------------------- if (Model%do_deep) then +!## CCPP ## GFS_DCNV_generic.F90/GFS_DCNV_generic_pre_run Note: The conditional +! above is not checked within the scheme, so the execution of the code below +! is controlled via its presence in the CCPP SDF. if (Model%do_ca) then do k=1,levs do i=1,im @@ -3519,10 +3684,11 @@ subroutine GFS_physics_driver & enddo enddo endif - +!*## CCPP ## if (.not. Model%ras .and. .not. Model%cscnv) then if (Model%imfdeepcnv == 1) then ! no random cloud top +!## CCPP ##* sascnvn.F/sascnvn_run call sascnvn (im, ix, levs, Model%jcap, dtp, del, & Statein%prsl, Statein%pgr, Statein%phil, clw(:,:,1:2), & Stateout%gq0(:,:,1), Stateout%gt0, Stateout%gu0, & @@ -3536,12 +3702,16 @@ subroutine GFS_physics_driver & Model%c1_deep, Model%betal_deep, Model%betas_deep, & Model%evfact_deep, Model%evfactl_deep, & Model%pgcon_deep) +!*## CCPP ## elseif (Model%imfdeepcnv == 2) then +!## CCPP ##* GFS_typedefs.f90/interstitial_setup_tracers if(.not. Model%satmedmf .and. .not. Model%trans_trac) then nsamftrac = 0 else nsamftrac = tottracer endif +!*## CCPP ## +!## CCPP ##* samfdeepcnv.f/samfdeepcnv_run call samfdeepcnv(im, ix, levs, dtp, itc, Model%ntchm, ntk, nsamftrac, & del, Statein%prsl, Statein%pgr, Statein%phil, clw, & Stateout%gq0(:,:,1), Stateout%gt0, & @@ -3557,6 +3727,7 @@ subroutine GFS_physics_driver & Model%c1_deep, Model%betal_deep, Model%betas_deep, & Model%evfact_deep, Model%evfactl_deep, & Model%pgcon_deep, Model%asolfac_deep) +!*## CCPP ## ! if (lprnt) print *,' rain1=',rain1(ipr) !elseif (Model%imfdeepcnv == 3) then ! if (Model%me==0) then @@ -3569,6 +3740,7 @@ subroutine GFS_physics_driver & ! stop ! end if elseif (Model%imfdeepcnv == 0) then ! random cloud top +!## CCPP ##* This is not in the CCPP yet. call sascnv (im, ix, levs, Model%jcap, dtp, del, & Statein%prsl, Statein%pgr, Statein%phil, clw(:,:,1:2), & Stateout%gq0(:,:,1), Stateout%gt0, Stateout%gu0, & @@ -3579,9 +3751,10 @@ subroutine GFS_physics_driver & ! QLCN, QICN, w_upi,cf_upi, CNV_MFD, CNV_PRC3, & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,imp_physics ) ! if (lprnt) print *,' rain1=',rain1(ipr),' rann=',rann(ipr,1) +!*## CCPP ## endif -! +!## CCPP ##* GFS_DCNV_generic.F90/GFS_DCNV_generic_post_run if (Model%npdf3d == 3 .and. Model%num_p3d == 4) then do k=1,levs do i=1,im @@ -3603,9 +3776,10 @@ subroutine GFS_physics_driver & if(Model%do_ca) then Coupling%cape(:) = cld1d(:) endif - +!*## CCPP ## ! else ! ras or cscnv +!## CCPP ## cs_conv.F90/cs_conv_pre_run fscav(:) = zero if (Model%cscnv) then ! Chikira-Sugiyama convection scheme (via CSU) @@ -3620,6 +3794,7 @@ subroutine GFS_physics_driver & ! dqdt(i,k,3) = clw(i,k,1) ! enddo ! enddo +!*## CCPP ## ! ! JLS NOTE: The convective mass fluxes (dt_mf, dd_mf and ud_mf) passed in and out of cs_conv have not been multiplied by @@ -3635,6 +3810,7 @@ subroutine GFS_physics_driver & ! JLS NOTE: The variable rain1 output from cs_convr (called prec inside the subroutine) is a precipitation flux (kg/m2/sec), ! not meters LWE like the other schemes. It is converted to m after the call to cs_convr. +!## CCPP ## cs_conv.F90/cs_conv_run call cs_convr (ix, im, levs, ntrac+1, nn, tottracer+3, & Model%nctp, otspt(1:ntrac+1,1:2), 1, & kdt, Stateout%gt0, Stateout%gq0(:,:,1:1), rain1, & @@ -3649,7 +3825,7 @@ subroutine GFS_physics_driver & w_upi, cf_upi, CNV_MFD, CNV_DQLDT, & ! w_upi, cf_upi, CNV_MFD, CNV_PRC3, CNV_DQLDT, & CLCN, CNV_FICE, CNV_NDROP, CNV_NICE, imp_physics) - +!*## CCPP ## ! if (lprnt) write(0,*)'aftcsgt0=',Stateout%gt0(ipr,:) ! if (lprnt) write(0,*)'aftcstke=',clw(ipr,1:25,ntk) @@ -3663,8 +3839,9 @@ subroutine GFS_physics_driver & ! &, Stateout%gq0(1:ix,1:levs,1),clw(1,1,2),clw(1,1,1) & ! &, ' cs_conv', grid%xlon(1:im), grid%xlat(1:im)) - +!## CCPP ##* Not in the CCPP. TODO: Does this need to be in cs_conv_post_run? rain1(:) = rain1(:) * (dtp*0.001) +!## CCPP ##* cs_conv.F90/cs_conv_post_run if (Model%do_aw) then do k=1,levs kk = min(k+1,levs) ! assuming no cloud top reaches the model top @@ -3673,7 +3850,7 @@ subroutine GFS_physics_driver & enddo enddo endif - +!*## CCPP ## ! if (lprnt) then ! write(0,*)' gt01=',stateout%gt0(ipr,:),' kdt=',kdt ! write(0,*)' gq01=',stateout%gq0(ipr,:,1),' kdt=',kdt @@ -3685,7 +3862,7 @@ subroutine GFS_physics_driver & else ! ras version 2 -! DH* this code not yet in CCPP (belongs to GFS_DCNV_generic_pre?) +!## CCPP ##* This code not yet in CCPP Note: Likely belongs in rascnv_pre. if (Model%ccwf(1) >= zero .or. Model%ccwf(2) >= 0) then do i=1,im ccwfac(i) = Model%ccwf(1)*work1(i) + Model%ccwf(2)*work2(i) @@ -3717,14 +3894,14 @@ subroutine GFS_physics_driver & ! if (ncld ==2) revap = .false. trcmin(:) = -999999.0 if (ntk-2 > 0) trcmin(ntk-2) = 1.0e-4 - +!*## CCPP ## ! if (lprnt) write(0,*)' gt04bras=',Stateout%gt0(ipr,:) ! if (lprnt) write(0,*)' gq04bras=',Stateout%gq0(ipr,:,1) ! if (lprnt) write(0,*)'befrasclw1=',clw(ipr,:,1) ! if (lprnt) write(0,*)'befrasclw2=',clw(ipr,:,2) ! if (lprnt) write(0,*)'befrastke=',clw(ipr,:,ntk) ! if (lprnt) write(0,*)'trcmin=',trcmin(ntk-2),' ntk=',ntk - +!## CCPP ## Not in CCPP yet. call rascnv (im, ix, levs, dtp, dtf, Tbd%rann, Stateout%gt0, & Stateout%gq0, Stateout%gu0, Stateout%gv0, clw, & tottracer, fscav, Statein%prsi, Statein%prsl, & @@ -3740,6 +3917,7 @@ subroutine GFS_physics_driver & CLCN, CNV_FICE, CNV_NDROP, CNV_NICE, imp_physics, & ! trcmin) trcmin, ntk) +!*## CCPP ## ! if (lprnt) write(0,*)' gt04=',Stateout%gt0(ipr,1:60) ! if (lprnt) write(0,*)' gq04=',Stateout%gq0(ipr,1:60,1) @@ -3759,10 +3937,13 @@ subroutine GFS_physics_driver & ! &,' cnv_prc3sum=',sum(cnv_prc3(ipr,1:levs)) ! if (lprnt) write(0,*)' gt04=',gt0(ipr,1:10) +!## CCPP ##* Not in CCPP yet. cld1d = 0 +!*## CCPP ## endif ! end if_not_ras +!## CCPP ##* GFS_DCNV_generic.F90/GFS_DCNV_generic_post if(Model%isppt_deep)then do k=1,levs do i=1,im @@ -3774,13 +3955,16 @@ subroutine GFS_physics_driver & enddo deallocate(savet, saveq, saveu, savev) endif - +!*## CCPP ## else ! no parameterized deep convection +!## CCPP ##* GFS_typedefs.F90/interstitial_phys_reset Note: These are only zeroed out +! initially, prior to calling physics. cld1d = zero rain1 = zero ud_mf = zero dd_mf = zero dt_mf = zero +!*## CCPP ## endif ! if (lprnt) then @@ -3790,6 +3974,7 @@ subroutine GFS_physics_driver & ! write(0,*)' aftcnvgq1=',(stateout%gq0(ipr,k,ntcw),k=1,levs) ! endif ! +!## CCPP ## GFS_DCNV_generic.F90/GFS_DCNV_generic_post_run do i=1,im Diag%rainc(i) = frain * rain1(i) enddo @@ -3815,7 +4000,9 @@ subroutine GFS_physics_driver & endif ! if (ldiag3d) endif ! end if_lssav +!*## CCPP ## ! +!## CCPP ##* This block not yet in CCPP. if (ldiag_ugwp) then tem = frain/dtp do k=1,levs @@ -3845,6 +4032,7 @@ subroutine GFS_physics_driver & ! from previous time step we need: LH-release + cld_top/bot + precip ! ! endif +!*## CCPP ## ! if (lprnt) write(7000,*)' bef cnvgwd gu0=',gu0(ipr,:) ! &,' lat=',lat,' kdt=',kdt,' me=',me @@ -3852,7 +4040,9 @@ subroutine GFS_physics_driver & ! !----------------Convective gravity wave drag parameterization starting -------- -! DH* this block is in gwdc_pre +!## CCPP ##* gwdc.f/gwdc_pre Note: The conditional above is not in the scheme, so +! the execution of the code below is controlled by its presence in the CCPP SDF +! --- ... calculate maximum convective heating rate if (Model%do_cnvgwd) then ! call convective gravity wave drag allocate(gwdcu(im,levs), gwdcv(im,levs)) @@ -3874,6 +4064,7 @@ subroutine GFS_physics_driver & do i=1,im if (work4(i) > zero) cumabs(i) = cumabs(i) / (dtp*work4(i)) enddo +!*## CCPP ## ! DH* 20180817 - note: the above non-CCPP code modifies work3, which until then was defined ! as the ratio of the exner function between midlayer and interface at lowest model layer: @@ -3952,6 +4143,9 @@ subroutine GFS_physics_driver & ! --- ... end check print ******************************************** +!## CCPP ##* gwdc.f/gwdc_run Note: The conditional above is not in the scheme, so +! the execution of the code below is controlled by its presence in the CCPP SDF + !GFDL replacing lat with "1" ! call gwdc(im, ix, im, levs, lat, gu0, gv0, gt0, gq0, dtp, & call gwdc (im, ix, im, levs, 1, Statein%ugrs, Statein%vgrs, & @@ -3982,7 +4176,7 @@ subroutine GFS_physics_driver & ! endif ! --- ... write out cloud top stress and wind tendencies - +!## CCPP ## gwdc.f/gwdc_post_run if (Model%lssav) then do i=1,im Diag%dugwd(i) = Diag%dugwd(i) + dusfcg(i)*dtf @@ -4013,6 +4207,7 @@ subroutine GFS_physics_driver & ! &gwdcu(ipr,k), ' gv0=', gv0(ipr,k),' gwdcv=',gwdcv(ipr,k) ! &,' k=',k enddo +!*## CCPP ## ! if (lprnt) then ! if (fhour >= fhourpr) then @@ -4042,6 +4237,7 @@ subroutine GFS_physics_driver & ! &,' lat=',lat,' kdt=',kdt,' me=',me !----------------Convective gravity wave drag parameterization over -------- +!## CCPP ## GFS_SCNV_generic.F90/GFS_SCNV_generic_pre_run if (Model%ldiag3d) then do k=1,levs do i=1,im @@ -4049,6 +4245,7 @@ subroutine GFS_physics_driver & enddo enddo endif +!*## CCPP ## if (.not. Model%do_shoc) then @@ -4056,6 +4253,7 @@ subroutine GFS_physics_driver & ! -------------------------------------- if (Model%imfshalcnv == 1) then ! opr option now at 2014 !----------------------- +!## CCPP ##* shalcnv.F/shalcnv_run call shalcnv (im, ix, levs, Model%jcap, dtp, del, Statein%prsl, & Statein%pgr, Statein%phil, clw, Stateout%gq0, & Stateout%gt0, Stateout%gu0, Stateout%gv0, rain1, & @@ -4063,7 +4261,9 @@ subroutine GFS_physics_driver & Diag%hpbl, hflx, evap, ud_mf, dt_mf, cnvw, cnvc, & Model%clam_shal, Model%c0s_shal, Model%c1_shal, & Model%pgcon_shal) +!*## CCPP ## +!## CCPP ##* GFS_SCNV_generic.F90/GFS_SCNV_generic_post_run do i=1,im Diag%rainc(i) = Diag%rainc(i) + frain * rain1(i) enddo @@ -4082,13 +4282,17 @@ subroutine GFS_physics_driver & enddo enddo endif +!*## CCPP ## elseif (Model%imfshalcnv == 2) then +!## CCPP ##* GFS_typedef.F90/interstitial_setup_tracers if(.not. Model%satmedmf .and. .not. Model%trans_trac) then nsamftrac = 0 else nsamftrac = tottracer endif +!*## CCPP ## +!## CCPP ##* samfshalcnv.f/samfshalcnv_run call samfshalcnv (im, ix, levs, dtp, itc, Model%ntchm, ntk, nsamftrac, & del, Statein%prsl, Statein%pgr, Statein%phil, clw, & Stateout%gq0(:,:,1), Stateout%gt0, & @@ -4098,7 +4302,8 @@ subroutine GFS_physics_driver & dt_mf, cnvw, cnvc, & Model%clam_shal, Model%c0s_shal, Model%c1_shal, & Model%pgcon_shal, Model%asolfac_shal) - +!*## CCPP ## +!## CCPP ##* GFS_SCNV_generic.F90/GFS_SCNV_generic_post_run do i=1,im Diag%rainc(i) = Diag%rainc(i) + frain * rain1(i) enddo @@ -4117,6 +4322,7 @@ subroutine GFS_physics_driver & enddo enddo endif +!*## CCPP ## !elseif (Model%imfshalcnv == 3) then !if (Model%me==0) write(0,*) "CCPP DEBUG: shallow convection of GF is called in gf_driver" @@ -4126,6 +4332,7 @@ subroutine GFS_physics_driver & elseif (Model%imfshalcnv == 0) then ! modified Tiedtke Shallow convecton !----------------------------------- +!## CCPP ## This block is not in the CCPP yet. levshc(:) = 0 do k=2,levs do i=1,im @@ -4154,8 +4361,10 @@ subroutine GFS_physics_driver & ! if (lprnt) print *,' levshcm=',levshcm,' gt0sha=',gt0(ipr,:) endif ! end if_imfshalcnv +!*## CCPP ## endif ! end if_shal_cnv +!## CCPP ## GFS_SCNV_generic.F90/GFS_SCNV_generic_post_run if (Model%lssav) then if (Model%ldiag3d) then do k=1,levs @@ -4181,7 +4390,7 @@ subroutine GFS_physics_driver & if (clw(i,k,2) <= -999.0) clw(i,k,2) = zero enddo enddo - +!*## CCPP ## ! if (lprnt) then ! write(0,*)' prsl=',prsl(ipr,:) ! write(0,*) ' del=',del(ipr,:) @@ -4189,6 +4398,7 @@ subroutine GFS_physics_driver & ! write(0,*) ' befshgq0=',gq0(ipr,:,1) ! endif +!## CCPP ##* gcm_shoc.F90/shoc_run elseif (Model%shocaftcnv) then ! if do_shoc is true and shocaftcnv is true call shoc if (imp_physics == Model%imp_physics_mg) then do k=1,levs @@ -4277,6 +4487,7 @@ subroutine GFS_physics_driver & ! write(0,*)' aft shoc gq0=',gq0(1,:,1),' lat=',lat ! write(0,*)' aft shoc gu0=',gu0(1,:),' lat=',lat ! +!*## CCPP ## endif ! if( .not. do_shoc) ! ! if (lprnt) then @@ -4286,6 +4497,7 @@ subroutine GFS_physics_driver & ! write(0,*) ' aftshgq0=',gq0(ipr,:,1) ! endif ! +!## CCPP ##* GFS_suite_interstitial.F90/GFS_suite_interstitial_4_run !------------------------------------------------------------------------------ ! --- update the tracers due to deep & shallow cumulus convective transport ! (except for suspended water and ice) @@ -4356,12 +4568,15 @@ subroutine GFS_physics_driver & enddo enddo endif ! end if_ntcw +!*## CCPP ## ! Legacy routine which determines convectve clouds - should be removed at some point - +!## CCPP ## cnvc90.f/cnvc90_run call cnvc90 (Model%clstp, im, ix, Diag%rainc, kbot, ktop, levs, Statein%prsi, & Tbd%acv, Tbd%acvb, Tbd%acvt, Cldprop%cv, Cldprop%cvb, Cldprop%cvt) +!*## CCPP ## +!## CCPP ##* This is not in the CCPP yet. if (Model%moist_adj) then ! moist convective adjustment ! --------------------------- ! @@ -4400,7 +4615,8 @@ subroutine GFS_physics_driver & ! endif ! endif endif ! moist convective adjustment over -! +!*## CCPP ## +!## CCPP ## GFS_MP_generic.F90/GFS_MP_generic_pre_run if (Model%ldiag3d .or. Model%do_aw) then do k=1,levs do i=1,im @@ -4412,8 +4628,10 @@ subroutine GFS_physics_driver & dqdt(1:im,:,n) = Stateout%gq0(1:im,:,n) enddo endif - +!*## CCPP ## ! dqdt_v : instaneous moisture tendency (kg/kg/sec) +!## CCPP ##* GFS_suite_interstitial.F90/GFS_suite_interstitial_4_run +! Note: ( these lines are relevant for shallow and deep convection) if (Model%cplchm) then do k=1,levs do i=1,im @@ -4421,15 +4639,16 @@ subroutine GFS_physics_driver & enddo enddo endif +!*## CCPP ## ! ! grid-scale condensation/precipitations and microphysics parameterization ! ------------------------------------------------------------------------ - +!## CCPP ##* This is not in the CCPP yet. if (ncld == 0) then ! no cloud microphysics call lrgscl (ix, im, levs, dtp, Stateout%gt0, Stateout%gq0, & Statein%prsl, del, Statein%prslk, rain1, clw) - +!*## CCPP ## else ! all microphysics if (imp_physics == Model%imp_physics_zhao_carr) then ! call zhao/carr/sundqvist microphysics ! ------------ @@ -4443,23 +4662,28 @@ subroutine GFS_physics_driver & ! write(0,*) ' beflsgw0=',gq0(ipr,:,3),' kdt=',kdt ! endif ! ------------------ +!## CCPP ##* This is not in the CCPP yet. if (Model%do_shoc) then call precpd_shoc (im, ix, levs, dtp, del, Statein%prsl, & Stateout%gq0(1,1,1), Stateout%gq0(1,1,ntcw), & Stateout%gt0, rain1, Diag%sr, rainp, rhc, & psautco_l, prautco_l, Model%evpco, Model%wminco, & Tbd%phy_f3d(1,1,ntot3d-2), lprnt, ipr) +!*## CCPP ## else +!## CCPP ##* gscond.f/zhaocarr_gscond_run call gscond (im, ix, levs, dtp, dtf, Statein%prsl, Statein%pgr, & Stateout%gq0(1,1,1), Stateout%gq0(1,1,ntcw), & Stateout%gt0, Tbd%phy_f3d(1,1,1), Tbd%phy_f3d(1,1,2), & Tbd%phy_f2d(1,1), Tbd%phy_f3d(1,1,3), & Tbd%phy_f3d(1,1,4), Tbd%phy_f2d(1,2), rhc,lprnt, ipr) - +!*## CCPP ## +!## CCPP ##* precpd.f/zhaocarr_precpd_run call precpd (im, ix, levs, dtp, del, Statein%prsl, & Stateout%gq0(1,1,1), Stateout%gq0(1,1,ntcw), & Stateout%gt0, rain1, Diag%sr, rainp, rhc, psautco_l, & prautco_l, Model%evpco, Model%wminco, lprnt, ipr) +!*## CCPP ## endif ! if (lprnt) then ! write(0,*)' prsl=',prsl(ipr,:) @@ -4472,6 +4696,7 @@ subroutine GFS_physics_driver & deallocate(rainp) elseif (imp_physics == Model%imp_physics_zhao_carr_pdf) then ! with pdf clouds +!## CCPP ##* These schemes are not in the CCPP yet. allocate(rainp(im,levs)) call gscondp (im, ix, levs, dtp, dtf, Statein%prsl, & Statein%pgr, Stateout%gq0(1,1,1), & @@ -4489,11 +4714,13 @@ subroutine GFS_physics_driver & Tbd%phy_f3d(1,1,Model%num_p3d+1), psautco_l, & prautco_l, Model%evpco, Model%wminco, lprnt, ipr) deallocate(rainp) +!*## CCPP ## ! if (lprnt) write(0,*) ' rain1=',rain1(ipr),' rainc=',rainc(ipr),' lat=',lat elseif (imp_physics == Model%imp_physics_thompson) then ! Thompson MP ! ------------ +!## CCPP ##* mp_thompson.F90/mp_thompson_run ims = 1 ; ime = ix ; kms = 1 ; kme = levs ; its = 1 ; ite = ix ; kts = 1 ; kte = levs if (Model%ltaerosol) then @@ -4532,9 +4759,9 @@ subroutine GFS_physics_driver & Diag%refl_10cm, Model%lradar, & Tbd%phy_f3d(:,:,1),Tbd%phy_f3d(:,:,2),Tbd%phy_f3d(:,:,3),me,Statein%phii) endif - +!*## CCPP ## elseif (imp_physics == Model%imp_physics_wsm6) then ! WSM6 - ! ----- +!## CCPP ##* This is not in the CCPP yet. ! ----- ims = 1 ; ime = ix ; kms = 1 ; kme = levs ; its = 1 ; ite = ix ; kts = 1 ; kte = levs call wsm6(Stateout%gt0, Statein%phii(1:im,1:levs+1), & @@ -4550,18 +4777,19 @@ subroutine GFS_physics_driver & Tbd%phy_f3d(:,:,1),Tbd%phy_f3d(:,:,2),Tbd%phy_f3d(:,:,3), & ims,ime, kms,kme, & its,ite, kts,kte) -! +!*## CCPP ## elseif (imp_physics == Model%imp_physics_mg) then ! MGB double-moment microphysics ! ------------------------------ - +!## CCPP ##* GFS_typedefs.F90/control_initialize kk = 5 if (Model%fprcp >= 2) kk = 6 - +!*## CCPP ## ! Acheng used clw here for other code to run smoothly and minimum change ! to make the code work. However, the nc and clw should be treated ! in other procceses too. August 28/2015; Hope that can be done next ! year. I believe this will make the physical interaction more reasonable ! Anning 12/5/2015 changed ntcw hold liquid only +!## CCPP ##* m_micro_insterstitial.F90/m_micro_pre_run if (Model%do_shoc) then skip_macro = Model%do_shoc if (Model%fprcp == 0) then @@ -4642,6 +4870,7 @@ subroutine GFS_physics_driver & Tbd%phy_f3d(i,k,1) = min(one, Tbd%phy_f3d(i,k,1) + clcn(i,k)) enddo enddo +!*## CCPP ## ! notice clw ix instead of im ! call m_micro_driver(im,ix,levs,flipv,del,dtp,prsl,prsi, @@ -4663,7 +4892,7 @@ subroutine GFS_physics_driver & ! do k=1,levs ! write(1000+me,*)' maxwatncb=',maxval(Stateout%gq0(1:im,k,ntlnc)),' k=',k,' kdt',kdt ! enddo - +!## CCPP ##* m_micro.F90/m_micro_run call m_micro_driver (im, ix, levs, Model%flipv, dtp, Statein%prsl, & Statein%prsi, Statein%phil, Statein%phii, & Statein%vvl, clw(1,1,2), QLCN, clw(1,1,1), QICN, & @@ -4688,6 +4917,7 @@ subroutine GFS_physics_driver & ! ipr, kdt, Grid%xlat, Grid%xlon) Model%mg_alf, Model%mg_qcmin, Model%pdfflag, & ipr, kdt, Grid%xlat, Grid%xlon, rhc) +!*## CCPP ## ! do k=1,levs ! write(1000+me,*)' maxwatnca=',maxval(Stateout%gq0(1:im,k,ntlnc)),' k=',k,' kdt=',kdt ! enddo @@ -4715,7 +4945,7 @@ subroutine GFS_physics_driver & ! if (lprnt) write(0,*)' qrna=',qrn(ipr,:),' kdt=',kdt ! if (lprnt) write(0,*)' qsnwa=',qsnw(ipr,:),' kdt=',kdt ! if (lprnt) write(0,*)' qglba',qgl(ipr,:),' kdt=',kdt - +!## CCPP ##* m_micro_interstitial.F90/m_micro_post_run tem = dtp * con_p001 / con_day if (abs(Model%fprcp) == 1 .or. mg3_as_mg2) then do k=1,levs @@ -4751,7 +4981,7 @@ subroutine GFS_physics_driver & Diag%snow(i) = tem * qsnw(i,1) Diag%graupel(i) = tem * qgl(i,1) enddo - +!*## CCPP ## endif ! if (lprnt) write(0,*)' cloudsm=',tbd%phy_f3d(ipr,:,1)*100,' kdt=',kdt @@ -4764,8 +4994,11 @@ subroutine GFS_physics_driver & elseif (imp_physics == Model%imp_physics_gfdl) then ! GFDL MP ! ------- do i = 1, im +!## CCPP ##* Not necessary in the CCPP. land (i,1) = frland(i) area (i,1) = Grid%area(i) +!*## CCPP ## +!## CCPP ##* gfdl_cloud_microphys.F90/gfdl_cloud_microphys_run rain0 (i,1) = zero snow0 (i,1) = zero ice0 (i,1) = zero @@ -4878,6 +5111,8 @@ subroutine GFS_physics_driver & enddo endif enddo +!*## CCPP ## +!## CCPP ##* maximum_hourly_diagnostics.F90/maximum_hourly_diagnsostics_run !Calculate hourly max 1-km agl and -10C reflectivity if (Model%lradar .and. & (imp_physics == Model%imp_physics_gfdl .or. & @@ -4898,7 +5133,8 @@ subroutine GFS_physics_driver & deallocate (refd) deallocate (refd263k) endif -! +!*## CCPP ## +!## CCPP ##* gfdl_cloud_microphys.F90/gfdl_cloud_microphys_run if(Model%effr_in) then call cloud_diagnosis (1, im, 1, levs, den(1:im,1:levs), & del(1:im,1:levs), islmsk(1:im), & @@ -4909,7 +5145,7 @@ subroutine GFS_physics_driver & Tbd%phy_f3d(1:im,1:levs,1), Tbd%phy_f3d(1:im,1:levs,2), & Tbd%phy_f3d(1:im,1:levs,3), Tbd%phy_f3d(1:im,1:levs,4), & Tbd%phy_f3d(1:im,1:levs,5)) - +!*## CCPP ## ! do k = 1, levs ! do i=1,im ! @@ -4944,6 +5180,9 @@ subroutine GFS_physics_driver & ! adjust sfc rainrate for conservation ! vertically integrate reduction of water increments, reduce precip by that amount +!## CCPP ##* cs_conv_aw_adj.F90/cs_conv_aw_adj_run Note: The conditional above +! is not checked in the scheme, so the control of the code below is through its +! inclusion in a CCPP SDF temrain1(:) = zero do k = 1,levs @@ -4982,9 +5221,34 @@ subroutine GFS_physics_driver & rain1(i) = max(rain1(i) - temrain1(i)*0.001, 0.0_kind_phys) enddo endif - +!*## CCPP ## +!## CCPP ##* GFS_MP_generic.F90/GFS_MP_generic_post_run Diag%rain(:) = Diag%rainc(:) + frain * rain1(:) - + +! --- get the amount of different precip type for Noah MP +! --- convert from m/dtp to mm/s + if (Model%lsm==Model%lsm_noahmp) then + if (Model%imp_physics == Model%imp_physics_mg .or. & + Model%imp_physics == Model%imp_physics_gfdl) then + !GJF: Should all precipitation rates have the same denominator below? + ! It appears that Diag%rain and Diag%rainc are on the dynamics time step, + ! but Diag%snow,graupel,ice are on the physics time step? This doesn't + ! matter as long as dtp=dtf (frain=1). + tem = 1.0 / (dtp*con_p001) + Sfcprop%draincprv(:) = tem * Diag%rainc(:) + Sfcprop%drainncprv(:) = tem * (frain * rain1(:)) + Sfcprop%dsnowprv(:) = tem * Diag%snow(:) + Sfcprop%dgraupelprv(:) = tem * Diag%graupel(:) + Sfcprop%diceprv(:) = tem * Diag%ice(:) + else + Sfcprop%draincprv(:) = 0.0 + Sfcprop%drainncprv(:) = 0.0 + Sfcprop%dsnowprv(:) = 0.0 + Sfcprop%dgraupelprv(:) = 0.0 + Sfcprop%diceprv(:) = 0.0 + endif + end if ! if (Model%lsm == Model%lsm_noahmp) + if (Model%cal_pre) then ! hchuang: add dominant precipitation type algorithm ! call calpreciptype (kdt, Model%nrcm, im, ix, levs, levs+1, & @@ -5051,7 +5315,8 @@ subroutine GFS_physics_driver & enddo endif endif - +!*## CCPP ## +!## CCPP ##* this block not yet in CCPP !-------------------------------- ! vay-2018 for Dycore-Tendencies save Stateout%X => Diag%dX3dt_cgw ! @@ -5061,6 +5326,8 @@ subroutine GFS_physics_driver & Diag%du3dt_cgw = Stateout%gu0 endif !-------------------------------- +!*## CCPP ## +!## CCPP ##* GFS_MP_generic.F90/GFS_MP_generic_post_run ! --- ... estimate t850 for rain-snow decision @@ -5149,9 +5416,10 @@ subroutine GFS_physics_driver & Coupling%rainc_cpl(i) = Coupling%rainc_cpl(i) + Diag%rainc(i) enddo endif - +!*## CCPP ## ! --- ... end coupling insertion +!## CCPP ##* GFS_surface_generic.F90/GFS_surface_generic_post_run ! --- ... total runoff is composed of drainage into water table and ! runoff at the surface and is accumulated in unit of meters if (Model%lssav) then @@ -5160,7 +5428,9 @@ subroutine GFS_physics_driver & Diag%srunoff(i) = Diag%srunoff(i) + runof(i) * dtf enddo endif - +!*## CCPP ## +!## CCPP ##* This block is not in the CCPP because data transfer between global +! and local variables is not necessary in the CCPP. ! --- ... return updated smsoil and stsoil to global arrays if (Model%frac_grid) then do k=1,lsoil @@ -5181,6 +5451,7 @@ subroutine GFS_physics_driver & enddo enddo endif +!*## CCPP ## ! --- ... calculate column precipitable water "pwat" Diag%pwat(:) = zero @@ -5237,7 +5508,8 @@ subroutine GFS_physics_driver & endif enddo endif - +!*## CCPP ## +!## CCPP ##* This block is not in the CCPP since it is not needed in the CCPP. deallocate (clw) if (allocated(cnvc)) deallocate(cnvc) if (allocated(cnvw)) deallocate(cnvw) @@ -5279,6 +5551,8 @@ subroutine GFS_physics_driver & w, qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt,p123,refl) deallocate (den) endif +!*## CCPP ## +!## CCPP ##* maximum_hourly_diagnostics.F90/maximum_hourly_diagnostics_run if (allocated(tke)) deallocate (tke) if (Model%cscnv) then deallocate (sigmatot, sigmafrac) @@ -5318,7 +5592,7 @@ subroutine GFS_physics_driver & Diag%T02MAX(I) = MAX(Diag%T02MAX(I), Sfcprop%t2m(i)) !<--- Hourly max 2m T Diag%T02MIN(I) = MIN(Diag%T02MIN(I), Sfcprop%t2m(i)) !<--- Hourly min 2m T enddo - +!*## CCPP ## ! if (kdt > 2 ) stop return !................................... diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index 69a7ebbb9..7a57f095e 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -159,10 +159,10 @@ subroutine FV3GFS_restart_write (IPD_Data, IPD_Restart, Atm_block, Model, fv_dom type(domain2d), intent(in) :: fv_domain character(len=32), optional, intent(in) :: timestamp - !--- read in surface data from chgres + !--- write surface data from chgres call sfc_prop_restart_write (IPD_Data%Sfcprop, Atm_block, Model, fv_domain, timestamp) - !--- read in physics restart data + !--- write physics restart data call phys_restart_write (IPD_Restart, Atm_block, Model, fv_domain, timestamp) end subroutine FV3GFS_restart_write @@ -287,7 +287,7 @@ subroutine FV3GFS_IPD_checksum (Model, IPD_Data, Atm_block) temp2d(i,j,55) = IPD_Data(nb)%Coupling%visbmui(ix) temp2d(i,j,56) = IPD_Data(nb)%Coupling%visdfui(ix) temp2d(i,j,57) = IPD_Data(nb)%Coupling%sfcdsw(ix) - temp2d(i,j,59) = IPD_Data(nb)%Coupling%sfcnsw(ix) + temp2d(i,j,58) = IPD_Data(nb)%Coupling%sfcnsw(ix) temp2d(i,j,59) = IPD_Data(nb)%Coupling%sfcdlw(ix) temp2d(i,j,60) = IPD_Data(nb)%Grid%xlon(ix) temp2d(i,j,61) = IPD_Data(nb)%Grid%xlat(ix) @@ -514,10 +514,18 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) nvar_s2o = 18 #ifdef CCPP if (Model%lsm == Model%lsm_ruc .and. warm_start) then - nvar_s2r = 6 + if(Model%rdlai) then + nvar_s2r = 7 + else + nvar_s2r = 6 + end if nvar_s3 = 5 else - nvar_s2r = 0 + if(Model%rdlai) then + nvar_s2r = 1 + else + nvar_s2r = 0 + endif nvar_s3 = 3 endif #else @@ -760,6 +768,11 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) sfc_name2(nvar_s2m+22) = 'tsnow' sfc_name2(nvar_s2m+23) = 'snowfall_acc' sfc_name2(nvar_s2m+24) = 'swe_snowfall_acc' + if (Model%rdlai) then + sfc_name2(nvar_s2m+25) = 'lai' + endif + else if (Model%lsm == Model%lsm_ruc .and. Model%rdlai) then + sfc_name2(nvar_s2m+19) = 'lai' #endif endif @@ -958,6 +971,11 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) Sfcprop(nb)%tsnow(ix) = sfc_var2(i,j,nvar_s2m+22) Sfcprop(nb)%snowfallac(ix) = sfc_var2(i,j,nvar_s2m+23) Sfcprop(nb)%acsnow(ix) = sfc_var2(i,j,nvar_s2m+24) + if (Model%rdlai) then + Sfcprop(nb)%xlaixy(ix) = sfc_var2(i,j,nvar_s2m+25) + endif + else if (Model%lsm == Model%lsm_ruc .and. Model%rdlai) then + Sfcprop(nb)%xlaixy(ix) = sfc_var2(i,j,nvar_s2m+19) elseif (Model%lsm == Model%lsm_noahmp) then !--- Extra Noah MP variables #else @@ -1085,8 +1103,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) do nb = 1, Atm_block%nblks do ix = 1, Atm_block%blksz(nb) Sfcprop(nb)%sncovr(ix) = 0.0 - !if (Sfcprop(nb)%slmsk(ix) > 0.001) then - if (Sfcprop(nb)%landfrac(ix) >= drythresh) then + if (Sfcprop(nb)%landfrac(ix) >= drythresh .or. Sfcprop(nb)%fice(ix) >= Model%min_seaice) then vegtyp = Sfcprop(nb)%vtype(ix) if (vegtyp == 0) vegtyp = 7 rsnow = 0.001*Sfcprop(nb)%weasd(ix)/snupx(vegtyp) @@ -1122,6 +1139,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) !#endif + if(Model%frac_grid) then ! 3-way composite do nb = 1, Atm_block%nblks do ix = 1, Atm_block%blksz(nb) if (Sfcprop(nb)%lakefrac(ix) > 0.0) then @@ -1135,8 +1153,21 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) if (Sfcprop(nb)%fice(ix) > 0. .and. Sfcprop(nb)%landfrac(ix)==0.) Sfcprop(nb)%slmsk(ix) = 2 ! land dominates over ice if co-exist enddo enddo + else !frac_grid=F + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + if (Sfcprop(nb)%lakefrac(ix) > 0.0) then + Sfcprop(nb)%oceanfrac(ix) = 0.0 ! lake & ocean don't coexist in a cell + if (Sfcprop(nb)%fice(ix) < Model%min_lakeice) Sfcprop(nb)%fice(ix) = 0. + else + Sfcprop(nb)%oceanfrac(ix) = 1.0 - Sfcprop(nb)%slmsk(ix) + if (Sfcprop(nb)%fice(ix) < Model%min_seaice) Sfcprop(nb)%fice(ix) = 0. + endif + enddo + enddo + end if - if(Model%frac_grid) then ! 3-way composite + if(Model%frac_grid) then ! 3-way composite do nb = 1, Atm_block%nblks do ix = 1, Atm_block%blksz(nb) Sfcprop(nb)%tsfco(ix) = max(con_tice, Sfcprop(nb)%tsfco(ix)) @@ -1211,7 +1242,6 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) Sfcprop(nb)%smoiseq(ix, 1:4) = missing_value Sfcprop(nb)%zsnsoxy(ix, -2:4) = missing_value - ! if (Sfcprop(nb)%slmsk(ix) > 0.01) then if (Sfcprop(nb)%landfrac(ix) >= drythresh) then Sfcprop(nb)%tvxy(ix) = Sfcprop(nb)%tsfcl(ix) @@ -1477,7 +1507,11 @@ subroutine sfc_prop_restart_write (Sfcprop, Atm_block, Model, fv_domain, timesta nvar2o = 18 #ifdef CCPP if (Model%lsm == Model%lsm_ruc) then - nvar2r = 6 + if (Model%rdlai) then + nvar2r = 7 + else + nvar2r = 6 + endif nvar3 = 5 else nvar2r = 0 @@ -1504,7 +1538,10 @@ subroutine sfc_prop_restart_write (Sfcprop, Atm_block, Model, fv_domain, timesta #ifdef CCPP if (Model%lsm == Model%lsm_ruc) then if (allocated(sfc_name2)) then - if (size(sfc_var3,dim=3).ne.Model%lsoil_lsm) then + ! Re-allocate if one or more of the dimensions don't match + if (size(sfc_name2).ne.nvar2m+nvar2o+nvar2mp+nvar2r .or. & + size(sfc_name3).ne.nvar3+nvar3mp .or. & + size(sfc_var3,dim=3).ne.Model%lsoil_lsm) then !--- deallocate containers and free restart container deallocate(sfc_name2) deallocate(sfc_name3) @@ -1611,6 +1648,9 @@ subroutine sfc_prop_restart_write (Sfcprop, Atm_block, Model, fv_domain, timesta sfc_name2(nvar2m+22) = 'tsnow' sfc_name2(nvar2m+23) = 'snowfall_acc' sfc_name2(nvar2m+24) = 'swe_snowfall_acc' + if (Model%rdlai) then + sfc_name2(nvar2m+25) = 'lai' + endif else if(Model%lsm == Model%lsm_noahmp) then #else ! Only needed when Noah MP LSM is used - 29 2D @@ -1684,17 +1724,19 @@ subroutine sfc_prop_restart_write (Sfcprop, Atm_block, Model, fv_domain, timesta #ifdef CCPP if (Model%lsm == Model%lsm_noah .or. Model%lsm == Model%lsm_noahmp) then - !--- names of the 2D variables to save + !--- names of the 3D variables to save sfc_name3(1) = 'stc' sfc_name3(2) = 'smc' sfc_name3(3) = 'slc' - sfc_name3(4) = 'snicexy' - sfc_name3(5) = 'snliqxy' - sfc_name3(6) = 'tsnoxy' - sfc_name3(7) = 'smoiseq' - sfc_name3(8) = 'zsnsoxy' + if (Model%lsm == Model%lsm_noahmp) then + sfc_name3(4) = 'snicexy' + sfc_name3(5) = 'snliqxy' + sfc_name3(6) = 'tsnoxy' + sfc_name3(7) = 'smoiseq' + sfc_name3(8) = 'zsnsoxy' + endif else if (Model%lsm == Model%lsm_ruc) then - !--- names of the 2D variables to save + !--- names of the 3D variables to save sfc_name3(1) = 'tslb' sfc_name3(2) = 'smois' sfc_name3(3) = 'sh2o' @@ -1813,6 +1855,9 @@ subroutine sfc_prop_restart_write (Sfcprop, Atm_block, Model, fv_domain, timesta sfc_var2(i,j,nvar2m+22) = Sfcprop(nb)%tsnow(ix) sfc_var2(i,j,nvar2m+23) = Sfcprop(nb)%snowfallac(ix) sfc_var2(i,j,nvar2m+24) = Sfcprop(nb)%acsnow(ix) + if (Model%rdlai) then + sfc_var2(i,j,nvar2m+25) = Sfcprop(nb)%xlaixy(ix) + endif else if (Model%lsm == Model%lsm_noahmp) then #else From 42bc16ea43b4e1602e4c1818b89cba404c486d81 Mon Sep 17 00:00:00 2001 From: "Shan.Sun" Date: Fri, 14 Feb 2020 16:38:37 +0000 Subject: [PATCH 06/15] Update .gitmodules and submodule pointer for ccpp-physics for code review and testing --- .gitmodules | 4 ++-- ccpp/physics | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.gitmodules b/.gitmodules index d253f6966..2689069c7 100644 --- a/.gitmodules +++ b/.gitmodules @@ -8,5 +8,5 @@ branch = master [submodule "ccpp/physics"] path = ccpp/physics - url = https://github.com/NCAR/ccpp-physics - branch = master + url = https://github.com/shansun6/ccpp-physics + branch = fractional_landmask diff --git a/ccpp/physics b/ccpp/physics index 01ed01fb0..f562f446e 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 01ed01fb0b3112e96eb619e0339d88fb0201982f +Subproject commit f562f446ed6c7bbc567c02df4c18fd98b1eb35b2 From e83422c6c04fe1d9e45e40a0c39d3f56ffee0c14 Mon Sep 17 00:00:00 2001 From: "Shan.Sun" Date: Fri, 14 Feb 2020 16:43:08 +0000 Subject: [PATCH 07/15] Revert change to smoiseq definition --- io/FV3GFS_io.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index 7a57f095e..6234ecdb7 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -1445,7 +1445,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) smc = smc - dx if ( abs (dx) < 1.e-6) exit enddo ! iteration - Sfcprop(nb)%smoiseq(ix,ns) = min(max(smc,drythresh),smcmax*0.99) + Sfcprop(nb)%smoiseq(ix,ns) = min(max(smc,1.e-4),smcmax*0.99) enddo ! ddz soil layer else ! bexp <= 0.0 Sfcprop(nb)%smoiseq(ix,1:4) = smcmax From b4eddd58e518ee5d25c873e059a15ead616d7213 Mon Sep 17 00:00:00 2001 From: "Shan.Sun" Date: Sat, 15 Feb 2020 18:57:39 +0000 Subject: [PATCH 08/15] Bug fix in FV3GFS_io.F90. Results are bitwise identical in standalone FV3 as well as coupled mode when frac_grid=F --- gfsphysics/GFS_layer/GFS_physics_driver.F90 | 5 +- io/FV3GFS_io.F90 | 87 +++++++++------------ 2 files changed, 38 insertions(+), 54 deletions(-) diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index e6e0f4c15..c4fdde1ef 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -1063,9 +1063,8 @@ subroutine GFS_physics_driver & !*## CCPP ## ! --- ... xw: transfer ice thickness & concentration from global to local variables !## CCPP ## global to local variable transfer not necessary for these two - Sfcprop%fice(i) = min(one,Sfcprop%fice(i)) ! sea/lake ice fraction wrt water portion, not whole cell - fice(i) = Sfcprop%fice(i) zice(i) = Sfcprop%hice(i) + fice(i) = Sfcprop%fice(i) !*## CCPP ##* !## CCPP ##* GFS_surface_composites.F90/GFS_surface_composites_pre_run tice(i) = Sfcprop%tisfc(i) @@ -1214,7 +1213,7 @@ subroutine GFS_physics_driver & enddo endif do i=1,im - if(wet(i) .or. icy(i)) then ! open water or ice + if(wet(i)) then ! Water zorl3(i,3) = Sfcprop%zorlo(i) tsfc3(i,3) = Sfcprop%tsfco(i) tsurf3(i,3) = Sfcprop%tsfco(i) diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index 6234ecdb7..f94a5eab9 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -935,6 +935,25 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) Sfcprop(nb)%tsfcl(ix) = sfc_var2(i,j,33) !--- sfcl (temp on land portion of a cell) Sfcprop(nb)%zorll(ix) = sfc_var2(i,j,34) !--- zorll (zorl on land portion of a cell) end if + + if(Model%frac_grid) then ! obtain landfrac from slmsk + Sfcprop(nb)%slmsk(ix) = ceiling(Sfcprop(nb)%landfrac(ix)) + if (Sfcprop(nb)%fice(ix) > 0. .and. Sfcprop(nb)%landfrac(ix)==0.) Sfcprop(nb)%slmsk(ix) = 2 ! land dominates ice if co-exist + else ! obtain slmsk from landfrac + if (Sfcprop(nb)%slmsk(ix) > 1.9) then + Sfcprop(nb)%landfrac(ix) = 0.0 + else + Sfcprop(nb)%landfrac(ix) = Sfcprop(nb)%slmsk(ix) + endif + end if + + if (Sfcprop(nb)%lakefrac(ix) > 0.0) then + Sfcprop(nb)%oceanfrac(ix) = 0.0 ! lake & ocean don't coexist in a cell + if (Sfcprop(nb)%fice(ix) < Model%min_lakeice) Sfcprop(nb)%fice(ix) = 0. + else + Sfcprop(nb)%oceanfrac(ix) = 1.0 - Sfcprop(nb)%landfrac(ix) + if (Sfcprop(nb)%fice(ix) < Model%min_seaice) Sfcprop(nb)%fice(ix) = 0. + endif ! !--- NSSTM variables if ((Model%nstf_name(1) > 0) .and. (Model%nstf_name(2) == 1)) then @@ -1117,57 +1136,29 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) enddo endif - if(Model%cplflx .or. Model%frac_grid) then - if (nint(sfc_var2(1,1,33)) == -9999) then - if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing tsfcl') - do nb = 1, Atm_block%nblks - do ix = 1, Atm_block%blksz(nb) - Sfcprop(nb)%tsfcl(ix) = Sfcprop(nb)%tsfco(ix) !--- compute tsfcl from existing variables + if(Model%cplflx .or. Model%frac_grid) then + if (nint(sfc_var2(1,1,33)) == -9999) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing tsfcl') + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + Sfcprop(nb)%tsfcl(ix) = Sfcprop(nb)%tsfco(ix) !--- compute tsfcl from existing variables + enddo enddo - enddo - endif + endif - if (nint(sfc_var2(1,1,34)) == -9999) then - if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorll') - do nb = 1, Atm_block%nblks - do ix = 1, Atm_block%blksz(nb) - Sfcprop(nb)%zorll(ix) = Sfcprop(nb)%zorlo(ix) !--- compute zorll from existing variables + if (nint(sfc_var2(1,1,34)) == -9999) then + if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorll') + do nb = 1, Atm_block%nblks + do ix = 1, Atm_block%blksz(nb) + Sfcprop(nb)%zorll(ix) = Sfcprop(nb)%zorlo(ix) !--- compute zorll from existing variables + enddo enddo - enddo + endif endif - endif !#endif - if(Model%frac_grid) then ! 3-way composite - do nb = 1, Atm_block%nblks - do ix = 1, Atm_block%blksz(nb) - if (Sfcprop(nb)%lakefrac(ix) > 0.0) then - Sfcprop(nb)%oceanfrac(ix) = 0.0 ! lake & ocean don't coexist in a cell - if (Sfcprop(nb)%fice(ix) < Model%min_lakeice) Sfcprop(nb)%fice(ix) = 0. - else - Sfcprop(nb)%oceanfrac(ix) = 1.0 - Sfcprop(nb)%landfrac(ix) !LHS:ocean frac [0:1] - if (Sfcprop(nb)%fice(ix) < Model%min_seaice) Sfcprop(nb)%fice(ix) = 0. - endif - Sfcprop(nb)%slmsk(ix) = ceiling(Sfcprop(nb)%landfrac(ix)) - if (Sfcprop(nb)%fice(ix) > 0. .and. Sfcprop(nb)%landfrac(ix)==0.) Sfcprop(nb)%slmsk(ix) = 2 ! land dominates over ice if co-exist - enddo - enddo - else !frac_grid=F - do nb = 1, Atm_block%nblks - do ix = 1, Atm_block%blksz(nb) - if (Sfcprop(nb)%lakefrac(ix) > 0.0) then - Sfcprop(nb)%oceanfrac(ix) = 0.0 ! lake & ocean don't coexist in a cell - if (Sfcprop(nb)%fice(ix) < Model%min_lakeice) Sfcprop(nb)%fice(ix) = 0. - else - Sfcprop(nb)%oceanfrac(ix) = 1.0 - Sfcprop(nb)%slmsk(ix) - if (Sfcprop(nb)%fice(ix) < Model%min_seaice) Sfcprop(nb)%fice(ix) = 0. - endif - enddo - enddo - end if - - if(Model%frac_grid) then ! 3-way composite + if(Model%frac_grid) then ! 3-way composite do nb = 1, Atm_block%nblks do ix = 1, Atm_block%blksz(nb) Sfcprop(nb)%tsfco(ix) = max(con_tice, Sfcprop(nb)%tsfco(ix)) @@ -1180,7 +1171,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) + Sfcprop(nb)%tsfco(ix) * (1.-Sfcprop(nb)%landfrac(ix)-tem) enddo enddo - else ! in this case ice fracion is fraction of water fraction + else do nb = 1, Atm_block%nblks do ix = 1, Atm_block%blksz(nb) !--- specify tsfcl/zorll from existing variable tsfco/zorlo @@ -1188,12 +1179,6 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) Sfcprop(nb)%zorll(ix) = Sfcprop(nb)%zorlo(ix) Sfcprop(nb)%zorl(ix) = Sfcprop(nb)%zorlo(ix) Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tsfco(ix) - Sfcprop(nb)%tsfco(ix) = max(con_tice, Sfcprop(nb)%tsfco(ix)) - if (Sfcprop(nb)%slmsk(ix) > 1.9) then - Sfcprop(nb)%landfrac(ix) = 0.0 - else - Sfcprop(nb)%landfrac(ix) = Sfcprop(nb)%slmsk(ix) - endif enddo enddo endif ! if (Model%frac_grid) From c13ed8710769ea530bc56054637115062bca2be3 Mon Sep 17 00:00:00 2001 From: "Shan.Sun" Date: Mon, 17 Feb 2020 21:14:34 +0000 Subject: [PATCH 09/15] Correcting comments --- io/FV3GFS_io.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index f94a5eab9..cf2589755 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -936,10 +936,10 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) Sfcprop(nb)%zorll(ix) = sfc_var2(i,j,34) !--- zorll (zorl on land portion of a cell) end if - if(Model%frac_grid) then ! obtain landfrac from slmsk + if(Model%frac_grid) then ! obtain slmsk from landfrac and fice Sfcprop(nb)%slmsk(ix) = ceiling(Sfcprop(nb)%landfrac(ix)) if (Sfcprop(nb)%fice(ix) > 0. .and. Sfcprop(nb)%landfrac(ix)==0.) Sfcprop(nb)%slmsk(ix) = 2 ! land dominates ice if co-exist - else ! obtain slmsk from landfrac + else ! obtain landfrac from slmsk if (Sfcprop(nb)%slmsk(ix) > 1.9) then Sfcprop(nb)%landfrac(ix) = 0.0 else From 22f0dcbf6b6d75b33c09d311123026b2dc2df63a Mon Sep 17 00:00:00 2001 From: "Shan.Sun" Date: Thu, 20 Feb 2020 22:09:03 +0000 Subject: [PATCH 10/15] Updating Coupling%tisfcin_cpl over ocean points only --- atmos_model.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/atmos_model.F90 b/atmos_model.F90 index e5d39c336..9bd70feb1 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -1674,7 +1674,9 @@ subroutine assign_importdata(rc) do i=isc,iec nb = Atm_block%blkno(i,j) ix = Atm_block%ixp(i,j) - IPD_Data(nb)%Coupling%tisfcin_cpl(ix) = datar8(i,j) + if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then + IPD_Data(nb)%Coupling%tisfcin_cpl(ix) = datar8(i,j) + endif enddo enddo endif From 13b04a49b9c08507b3917510e96f3774d25edaaa Mon Sep 17 00:00:00 2001 From: "Shan.Sun" Date: Thu, 12 Mar 2020 17:31:33 +0000 Subject: [PATCH 11/15] Merging develop (fc9a4c2) into branch fractional_landmask --- atmos_model.F90 | 96 +++- gfsphysics/GFS_layer/GFS_physics_driver.F90 | 536 +++++++++++--------- io/FV3GFS_io.F90 | 10 +- 3 files changed, 382 insertions(+), 260 deletions(-) diff --git a/atmos_model.F90 b/atmos_model.F90 index 9bd70feb1..b9a91b94c 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -223,7 +223,8 @@ module atmos_model_mod #endif real(kind=IPD_kind_phys), parameter :: zero = 0.0_IPD_kind_phys, & - one = 1.0_IPD_kind_phys + one = 1.0_IPD_kind_phys, & + epsln = 1.0e-10_IPD_kind_phys contains @@ -1596,6 +1597,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 + real(kind=IPD_kind_phys) :: tem logical found, isFieldCreated, lcpl_fice ! !------------------------------------------------------------------------------ @@ -1663,6 +1665,29 @@ subroutine assign_importdata(rc) ! endif ! endif + +! get sea-state dependent surface roughness (if cplwav2atm=true) +!---------------------------- + fldname = 'wave_z0_roughness_length' + if (trim(impfield_name) == trim(fldname)) then + findex = QueryFieldList(ImportFieldsList,fldname) + if (importFieldsValid(findex) .and. IPD_control%cplwav2atm) then +!$omp parallel do default(shared) private(i,j,nb,ix) + do j=jsc,jec + do i=isc,iec + nb = Atm_block%blkno(i,j) + ix = Atm_block%ixp(i,j) + if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then + tem = 100.0 * max(zero, min(0.1, datar8(i,j))) + IPD_Data(nb)%Coupling%zorlwav_cpl(ix) = tem + IPD_Data(nb)%Sfcprop%zorlo(ix) = tem + + endif + enddo + enddo + endif + endif + ! get sea ice surface temperature !-------------------------------- fldname = 'sea_ice_surface_temperature' @@ -1699,7 +1724,6 @@ subroutine assign_importdata(rc) if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then IPD_Data(nb)%Coupling%tseain_cpl(ix) = datar8(i,j) IPD_Data(nb)%Sfcprop%tsfco(ix) = datar8(i,j) -! IPD_Data(nb)%Sfcprop%tsfc(ix) = datar8(i,j) endif enddo enddo @@ -1720,15 +1744,18 @@ subroutine assign_importdata(rc) nb = Atm_block%blkno(i,j) ix = Atm_block%ixp(i,j) IPD_Data(nb)%Coupling%ficein_cpl(ix) = zero + IPD_Data(nb)%Coupling%slimskin_cpl(ix) = IPD_Data(nb)%Sfcprop%slmsk(ix) if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then IPD_Data(nb)%Coupling%ficein_cpl(ix) = max(zero, min(one, datar8(i,j)/IPD_Data(nb)%Sfcprop%oceanfrac(ix))) !LHS: ice frac wrt water area if (IPD_Data(nb)%Coupling%ficein_cpl(ix) >= IPD_control%min_seaice) then - if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) == one) IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points + if (abs(one-IPD_Data(nb)%Sfcprop%oceanfrac(ix)) < epsln) IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points IPD_Data(nb)%Coupling%slimskin_cpl(ix) = 4. else - if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) == one) IPD_Data(nb)%Sfcprop%slmsk(ix) = zero - IPD_Data(nb)%Coupling%slimskin_cpl(ix) = zero IPD_Data(nb)%Coupling%ficein_cpl(ix) = zero + if (abs(one-IPD_Data(nb)%Sfcprop%oceanfrac(ix)) < epsln) then + IPD_Data(nb)%Sfcprop%slmsk(ix) = zero + IPD_Data(nb)%Coupling%slimskin_cpl(ix) = zero + end if endif endif enddo @@ -1739,7 +1766,7 @@ subroutine assign_importdata(rc) ! get upward LW flux: for sea ice covered area !---------------------------------------------- - fldname = 'mean_up_lw_flx' + fldname = 'mean_up_lw_flx_ice' if (trim(impfield_name) == trim(fldname)) then findex = QueryFieldList(ImportFieldsList,fldname) if (importFieldsValid(findex)) then @@ -1766,7 +1793,7 @@ subroutine assign_importdata(rc) ! get latent heat flux: for sea ice covered area !------------------------------------------------ - fldname = 'mean_laten_heat_flx' + fldname = 'mean_laten_heat_flx_atm_into_ice' if (trim(impfield_name) == trim(fldname)) then findex = QueryFieldList(ImportFieldsList,fldname) if (importFieldsValid(findex)) then @@ -1786,7 +1813,7 @@ subroutine assign_importdata(rc) ! get sensible heat flux: for sea ice covered area !-------------------------------------------------- - fldname = 'mean_sensi_heat_flx' + fldname = 'mean_sensi_heat_flx_atm_into_ice' if (trim(impfield_name) == trim(fldname)) then findex = QueryFieldList(ImportFieldsList,fldname) if (importFieldsValid(findex)) then @@ -1806,7 +1833,7 @@ subroutine assign_importdata(rc) ! get zonal compt of momentum flux: for sea ice covered area !------------------------------------------------------------ - fldname = 'mean_zonal_moment_flx' + fldname = 'stress_on_air_ice_zonal' if (trim(impfield_name) == trim(fldname)) then findex = QueryFieldList(ImportFieldsList,fldname) if (importFieldsValid(findex)) then @@ -1826,7 +1853,7 @@ subroutine assign_importdata(rc) ! get meridional compt of momentum flux: for sea ice covered area !----------------------------------------------------------------- - fldname = 'mean_merid_moment_flx' + fldname = 'stress_on_air_ice_merid' if (trim(impfield_name) == trim(fldname)) then findex = QueryFieldList(ImportFieldsList,fldname) if (importFieldsValid(findex)) then @@ -1905,6 +1932,7 @@ subroutine assign_importdata(rc) IPD_Data(nb)%Sfcprop%hice(ix) = IPD_Data(nb)%Coupling%hicein_cpl(ix) IPD_Data(nb)%Sfcprop%snowd(ix) = IPD_Data(nb)%Coupling%hsnoin_cpl(ix) else + IPD_Data(nb)%Sfcprop%tisfc(ix) = IPD_Data(nb)%Coupling%tseain_cpl(ix) IPD_Data(nb)%Sfcprop%fice(ix) = zero IPD_Data(nb)%Sfcprop%hice(ix) = zero IPD_Data(nb)%Sfcprop%snowd(ix) = zero @@ -1915,12 +1943,30 @@ subroutine assign_importdata(rc) IPD_Data(nb)%Coupling%dvsfcin_cpl(ix) = -99999.0 ! ,, IPD_Data(nb)%Coupling%dtsfcin_cpl(ix) = -99999.0 ! ,, IPD_Data(nb)%Coupling%ulwsfcin_cpl(ix) = -99999.0 ! ,, - if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) == one) IPD_Data(nb)%Sfcprop%slmsk(ix) = zero ! 100% open water + if (abs(one-IPD_Data(nb)%Sfcprop%oceanfrac(ix)) < epsln) & + IPD_Data(nb)%Coupling%slimskin_cpl(ix) = zero ! 100% open water endif endif enddo enddo endif +! +!------------------------------------------------------------------------------- +! do j=jsc,jec +! do i=isc,iec +! nb = Atm_block%blkno(i,j) +! ix = Atm_block%ixp(i,j) +! if (abs(IPD_Data(nb)%Grid%xlon_d(ix)-2.89) < 0.1 .and. & +! abs(IPD_Data(nb)%Grid%xlat_d(ix)+58.99) < 0.1) then +! write(0,*)' in assign tisfc=',IPD_Data(nb)%Sfcprop%tisfc(ix), & +! ' oceanfrac=',IPD_Data(nb)%Sfcprop%oceanfrac(ix),' i=',i,' j=',j,& +! ' tisfcin=',IPD_Data(nb)%Coupling%tisfcin_cpl(ix), & +! ' fice=',IPD_Data(nb)%Sfcprop%fice(ix) +! endif +! enddo +! enddo +!------------------------------------------------------------------------------- +! rc=0 ! @@ -2517,7 +2563,7 @@ subroutine setup_exportdata (rc) exportData(i,j,idx) = DYCORE_Data(nb)%coupling%t_bot(ix) else exportData(i,j,idx) = zero - endif + endif enddo enddo endif @@ -2535,7 +2581,7 @@ subroutine setup_exportdata (rc) exportData(i,j,idx) = DYCORE_Data(nb)%coupling%tr_bot(ix,1) else exportData(i,j,idx) = zero - endif + endif enddo enddo endif @@ -2586,7 +2632,7 @@ subroutine setup_exportdata (rc) exportData(i,j,idx) = DYCORE_Data(nb)%coupling%p_bot(ix) else exportData(i,j,idx) = zero - endif + endif enddo enddo endif @@ -2602,8 +2648,8 @@ subroutine setup_exportdata (rc) if (associated(DYCORE_Data(nb)%coupling%z_bot)) then exportData(i,j,idx) = DYCORE_Data(nb)%coupling%z_bot(ix) else - exportData(i,j,idx) = zero - endif + exportData(i,j,idx) = zero + endif enddo enddo endif @@ -2622,14 +2668,14 @@ subroutine setup_exportdata (rc) enddo enddo endif - endif !cplflx + endif !cplflx !--- ! Fill the export Fields for ESMF/NUOPC style coupling call fillExportFields(exportData) !--- - if (IPD_Control%cplflx) then + if (IPD_Control%cplflx) then ! zero out accumulated fields !$omp parallel do default(shared) private(i,j,nb,ix) do j=jsc,jec @@ -2662,12 +2708,12 @@ subroutine setup_exportdata (rc) end subroutine setup_exportdata - subroutine addLsmask2grid(fcstgrid, rc) + subroutine addLsmask2grid(fcstGrid, rc) use ESMF ! implicit none - type(ESMF_Grid) :: fcstgrid + type(ESMF_Grid) :: fcstGrid integer, optional, intent(out) :: rc ! ! local vars @@ -2675,7 +2721,7 @@ subroutine addLsmask2grid(fcstgrid, rc) integer i, j, nb, ix ! integer CLbnd(2), CUbnd(2), CCount(2), TLbnd(2), TUbnd(2), TCount(2) type(ESMF_StaggerLoc) :: staggerloc - integer, allocatable :: lsmask(:,:) + integer, allocatable :: lsmask(:,:) integer(kind=ESMF_KIND_I4), pointer :: maskPtr(:,:) ! isc = IPD_control%isc @@ -2690,16 +2736,16 @@ subroutine addLsmask2grid(fcstgrid, rc) nb = Atm_block%blkno(i,j) ix = Atm_block%ixp(i,j) ! use land sea mask: land:1, ocean:0 - lsmask(i,j) = floor(IPD_Data(nb)%SfcProp%landfrac(ix)) + lsmask(i,j) = floor(one + epsln - IPD_Data(nb)%SfcProp%oceanfrac(ix)) enddo enddo ! ! Get mask - call ESMF_GridAddItem(fcstgrid, itemflag=ESMF_GRIDITEM_MASK, & + call ESMF_GridAddItem(fcstGrid, itemflag=ESMF_GRIDITEM_MASK, & staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return -! call ESMF_GridGetItemBounds(fcstgrid, itemflag=ESMF_GRIDITEM_MASK, & +! call ESMF_GridGetItemBounds(fcstGrid, itemflag=ESMF_GRIDITEM_MASK, & ! staggerloc=ESMF_STAGGERLOC_CENTER, computationalLBound=ClBnd, & ! computationalUBound=CUbnd, computationalCount=Ccount, & ! totalLBound=TLbnd, totalUBound=TUbnd, totalCount=Tcount, rc=rc) @@ -2708,7 +2754,7 @@ subroutine addLsmask2grid(fcstgrid, rc) ! 'TlBnd=',TlBnd,'TUbnd=',TUbnd,'Tcount=',Tcount ! if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - call ESMF_GridGetItem(fcstgrid, itemflag=ESMF_GRIDITEM_MASK, & + call ESMF_GridGetItem(fcstGrid, itemflag=ESMF_GRIDITEM_MASK, & staggerloc=ESMF_STAGGERLOC_CENTER,farrayPtr=maskPtr, rc=rc) ! print *,'in set up grid, aft get maskptr, rc=',rc, 'size=',size(maskPtr,1),size(maskPtr,2), & ! 'bound(maskPtr)=', LBOUND(maskPtr,1),LBOUND(maskPtr,2),UBOUND(maskPtr,1),UBOUND(maskPtr,2) diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index c4fdde1ef..22bc0e832 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -43,6 +43,7 @@ module module_physics_driver !--- CONSTANT PARAMETERS real(kind=kind_phys), parameter :: hocp = con_hvap/con_cp + real(kind=kind_phys), parameter :: epsln = 1.0d-10 real(kind=kind_phys), parameter :: qmin = 1.0d-10 real(kind=kind_phys), parameter :: qsmall = 1.0d-20 real(kind=kind_phys), parameter :: rainmin = 1.0d-13 @@ -829,8 +830,12 @@ subroutine GFS_physics_driver & ! lprnt = .false. ! do i=1,im -! lprnt = kdt >= 1 .and. abs(grid%xlon(i)*rad2dg-29.55) < 0.201 & -! .and. abs(grid%xlat(i)*rad2dg+59.62) < 0.201 +! lprnt = kdt >= 1 .and. abs(grid%xlon(i)*rad2dg-97.50) < 0.101 & +! .and. abs(grid%xlat(i)*rad2dg-24.48) < 0.101 +! lprnt = kdt >= 1 .and. abs(grid%xlon(i)*rad2dg-293.91) < 0.101 & +! .and. abs(grid%xlat(i)*rad2dg+72.02) < 0.101 +! lprnt = kdt >= 1 .and. abs(grid%xlon(i)*rad2dg-113.48) < 0.101 & +! .and. abs(grid%xlat(i)*rad2dg-21.07) < 0.101 ! lprnt = kdt >= 1 .and. abs(grid%xlon(i)*rad2dg-169.453) < 0.501 & ! .and. abs(grid%xlat(i)*rad2dg-72.96) < 0.501 ! if (kdt == 1) & @@ -1139,9 +1144,10 @@ subroutine GFS_physics_driver & fice(i) = zero endif endif - if (fice(i) < one ) then - wet(i)=.true. !there is some open ocean/lake water! + if (fice(i) < one) then + wet(i)=.true. !some open ocean/lake water exists if (.not. Model%cplflx) Sfcprop%tsfco(i) = max(Sfcprop%tsfco(i), Sfcprop%tisfc(i), tgice) + end if else fice(i) = zero @@ -1206,11 +1212,16 @@ subroutine GFS_physics_driver & enddo if (.not. Model%cplflx .or. .not. Model%frac_grid) then - do i=1,im - Sfcprop%zorll(i) = Sfcprop%zorl(i) - Sfcprop%zorlo(i) = Sfcprop%zorl(i) -! Sfcprop%tisfc(i) = Sfcprop%tsfc(i) - enddo + if (Model%cplwav2atm) then + do i=1,im + Sfcprop%zorll(i) = Sfcprop%zorl(i) + enddo + else + do i=1,im + Sfcprop%zorll(i) = Sfcprop%zorl(i) + Sfcprop%zorlo(i) = Sfcprop%zorl(i) + enddo + endif endif do i=1,im if(wet(i)) then ! Water @@ -1705,10 +1716,10 @@ subroutine GFS_physics_driver & tsurf3(i,3) = tsurf3(i,3) + tem endif enddo - if (Model%cplflx) then + if (Model%cplflx) then ! apply only at ocean points tem1 = half / omz1 do i=1,im - if (wet(i)) then + if (wet(i) .and. Sfcprop%oceanfrac(i) > zero) then tem2 = one / Sfcprop%xz(i) dt_warm = (Sfcprop%xt(i)+Sfcprop%xt(i)) * tem2 if ( Sfcprop%xz(i) > omz1) then @@ -1719,7 +1730,7 @@ subroutine GFS_physics_driver & - Sfcprop%z_c(i)*Sfcprop%dt_cool(i))*tem1 endif TSEAl(i) = Sfcprop%tref(i) + dt_warm - Sfcprop%dt_cool(i) -! - (Sfcprop%oro(i)-Sfcprop%oro_uf(i))*rlapse +! - (Sfcprop%oro(i)-Sfcprop%oro_uf(i))*rlapse tsurf3(i,3) = TSEAl(i) endif enddo @@ -1766,13 +1777,12 @@ subroutine GFS_physics_driver & zsea1 = 0.001*real(Model%nstf_name(4)) zsea2 = 0.001*real(Model%nstf_name(5)) call get_dtzm_2d (Sfcprop%xt, Sfcprop%xz, Sfcprop%dt_cool, & - Sfcprop%z_c, wet, zsea1, zsea2, & - im, 1, dtzm) + Sfcprop%z_c, wet, zsea1, zsea2, im, 1, dtzm) do i=1,im ! if (wet(i) .and. .not.icy(i)) then ! if (wet(i) .and. (Model%frac_grid .or. .not. icy(i))) then if (wet(i)) then - tsfc3(i,3) = max(271.2,Sfcprop%tref(i) + dtzm(i)) + tsfc3(i,3) = max(tgice,Sfcprop%tref(i) + dtzm(i)) ! tsfc3(i,3) = max(271.2,Sfcprop%tref(i) + dtzm(i)) - & ! (Sfcprop%oro(i)-Sfcprop%oro_uf(i))*rlapse endif @@ -2028,8 +2038,8 @@ subroutine GFS_physics_driver & Sfcprop%tsfc(i) = txl*tsfc3(i,1) + txi*tice(i) + txo*tsfc3(i,3) ! Sfcprop%tsfc(i) = txl*tsfc3(i,1) + txi*tsfc3(i,2) + txo*tsfc3(i,3) - Diag%cmm(i) = txl*cmm3(i,1) + txi*cmm3(i,2) + txo*cmm3(i,3) - Diag%chh(i) = txl*chh3(i,1) + txi*chh3(i,2) + txo*chh3(i,3) +! Diag%cmm(i) = txl*cmm3(i,1) + txi*cmm3(i,2) + txo*cmm3(i,3) +! Diag%chh(i) = txl*chh3(i,1) + txi*chh3(i,2) + txo*chh3(i,3) Sfcprop%zorll(i) = zorl3(i,1) Sfcprop%zorlo(i) = zorl3(i,3) @@ -2105,7 +2115,7 @@ subroutine GFS_physics_driver & Sfcprop%zorll(i) = zorl3(i,1) Sfcprop%zorlo(i) = zorl3(i,3) - if (flag_cice(i)) then ! this was already done for lake ice in sfc_sice + if (flag_cice(i) .and. wet(i)) then ! this was already done for lake ice in sfc_sice txi = fice(i) txo = one - txi evap(i) = txi * evap3(i,2) + txo * evap3(i,3) @@ -2201,8 +2211,8 @@ subroutine GFS_physics_driver & Coupling%nlwsfc_cpl (i) = Coupling%nlwsfc_cpl(i) + Coupling%nlwsfci_cpl(i)*dtf Coupling%t2mi_cpl (i) = Sfcprop%t2m(i) Coupling%q2mi_cpl (i) = Sfcprop%q2m(i) -! Coupling%tsfci_cpl (i) = Sfcprop%tsfc(i) - Coupling%tsfci_cpl (i) = tsfc3(i,3) + Coupling%tsfci_cpl (i) = Sfcprop%tsfc(i) +! Coupling%tsfci_cpl (i) = tsfc3(i,3) Coupling%psurfi_cpl (i) = Statein%pgr(i) enddo @@ -2330,12 +2340,16 @@ subroutine GFS_physics_driver & dvsfc1, dtsfc1, dqsfc1, dkt, Diag%hpbl, kinver, & Model%xkzm_m, Model%xkzm_h, Model%xkzm_s, Model%xkzminv, & lprnt, ipr, me) -! if (lprnt) write(0,*)'aftmonshoc=',Statein%tgrs(ipr,:) -! if (lprnt) write(0,*)'aftmonshocq=',Statein%qgrs(ipr,:,1) -! if (lprnt) write(0,*)'aftmonshoctke=',Statein%qgrs(ipr,:,ntke) -! if (lprnt) write(0,*)'aftmonice=',Statein%qgrs(ipr,:,ntiw) -! if (lprnt) write(0,*)'aftmonwat=',Statein%qgrs(ipr,:,ntcw) -! if (lprnt) write(0,*)'aftmonshocdtdt=',dtdt(ipr,1:10) +! if (lprnt) then +! write(0,*)' aftpbl dtdt=',dtdt(ipr,:) +! write(0,*)' aftpbl dqdtv=',dqdt(ipr,:,1) +! write(0,*)'aftmonshoc=',Statein%tgrs(ipr,:) +! write(0,*)'aftmonshocq=',Statein%qgrs(ipr,:,1) +! write(0,*)'aftmonshoctke=',Statein%qgrs(ipr,:,ntke) +! write(0,*)'aftmonice=',Statein%qgrs(ipr,:,ntiw) +! write(0,*)'aftmonwat=',Statein%qgrs(ipr,:,ntcw) +! write(0,*)'aftmonshocdtdt=',dtdt(ipr,1:10) +! endif else if (Model%satmedmf) then if (Model%isatmedmf == 0) then ! initial version of satmedmfvdif (Nov 2018) @@ -2810,7 +2824,8 @@ subroutine GFS_physics_driver & endif ! if (lprnt) then -! write(0,*) ' dusfc1=',dusfc1(ipr),' kdt=',kdt,' lat=',lat +! write(0,*) ' dusfc1=',dusfc1(ipr),' kdt=',kdt +! write(0,*) ' dvsfc1=',dvsfc1(ipr),' kdt=',kdt ! write(0,*)' dtsfc1=',dtsfc1(ipr) ! write(0,*)' dqsfc1=',dqsfc1(ipr) ! write(0,*)' dtdtc=',(dtdt(ipr,k),k=1,15) @@ -2823,34 +2838,31 @@ subroutine GFS_physics_driver & if (Model%cplflx) then do i=1,im if (Sfcprop%oceanfrac(i) > zero) then ! Ocean only, NO LAKES -! if (Sfcprop%fice(i) == Sfcprop%oceanfrac(i)) then ! use results from CICE -! Coupling%dusfci_cpl(i) = dusfc_cice(i) -! Coupling%dvsfci_cpl(i) = dvsfc_cice(i) -! Coupling%dtsfci_cpl(i) = dtsfc_cice(i) -! Coupling%dqsfci_cpl(i) = dqsfc_cice(i) -! elseif (dry(i) .or. icy(i)) then ! use stress_ocean from sfc_diff for opw component at mixed point - if (wet(i)) then ! use stress_ocean from sfc_diff for opw component at mixed point - if (icy(i) .or. dry(i)) then - tem1 = max(Diag%q1(i), 1.e-8) - rho = Statein%prsl(i,1) / (con_rd*Diag%t1(i)*(one+con_fvirt*tem1)) - if (wind(i) > zero) then - tem = - rho * stress3(i,3) / wind(i) - Coupling%dusfci_cpl(i) = tem * Statein%ugrs(i,1) ! U-momentum flux - Coupling%dvsfci_cpl(i) = tem * Statein%vgrs(i,1) ! V-momentum flux - else - Coupling%dusfci_cpl(i) = zero - Coupling%dvsfci_cpl(i) = zero - endif - Coupling%dtsfci_cpl(i) = con_cp * rho * hflx3(i,3) ! sensible heat flux over open ocean - Coupling%dqsfci_cpl(i) = con_hvap * rho * evap3(i,3) ! latent heat flux over open ocean + if (Sfcprop%fice(i) >= one - epsln) then ! no open water, thus use results from CICE + Coupling%dusfci_cpl(i) = dusfc_cice(i) + Coupling%dvsfci_cpl(i) = dvsfc_cice(i) + Coupling%dtsfci_cpl(i) = dtsfc_cice(i) + Coupling%dqsfci_cpl(i) = dqsfc_cice(i) + elseif (icy(i) .or. dry(i)) then ! use stress_ocean from sfc_diff for opw component at mixed point + tem1 = max(Diag%q1(i), 1.e-8) + rho = Statein%prsl(i,1) / (con_rd*Diag%t1(i)*(one+con_fvirt*tem1)) + if (wind(i) > zero) then + tem = - rho * stress3(i,3) / wind(i) + Coupling%dusfci_cpl(i) = tem * Statein%ugrs(i,1) ! U-momentum flux + Coupling%dvsfci_cpl(i) = tem * Statein%vgrs(i,1) ! V-momentum flux + else + Coupling%dusfci_cpl(i) = zero + Coupling%dvsfci_cpl(i) = zero + endif + Coupling%dtsfci_cpl(i) = con_cp * rho * hflx3(i,3) ! sensible heat flux over open ocean + Coupling%dqsfci_cpl(i) = con_hvap * rho * evap3(i,3) ! latent heat flux over open ocean ! if (lprnt .and. i == ipr) write(0,*)' hflx33=',hflx3(i,3),' evap33=',evap3(i,3), & ! ' con_cp=',con_cp,' rho=',rho,' con_hvap=',con_hvap - else ! use results from PBL scheme for 100% open ocean - Coupling%dusfci_cpl(i) = dusfc1(i) - Coupling%dvsfci_cpl(i) = dvsfc1(i) - Coupling%dtsfci_cpl(i) = dtsfc1(i) - Coupling%dqsfci_cpl(i) = dqsfc1(i) - endif + else ! use results from PBL scheme for 100% open ocean + Coupling%dusfci_cpl(i) = dusfc1(i) + Coupling%dvsfci_cpl(i) = dvsfc1(i) + Coupling%dtsfci_cpl(i) = dtsfc1(i) + Coupling%dqsfci_cpl(i) = dqsfc1(i) endif Coupling%dusfc_cpl (i) = Coupling%dusfc_cpl(i) + Coupling%dusfci_cpl(i) * dtf @@ -2858,6 +2870,11 @@ subroutine GFS_physics_driver & Coupling%dtsfc_cpl (i) = Coupling%dtsfc_cpl(i) + Coupling%dtsfci_cpl(i) * dtf Coupling%dqsfc_cpl (i) = Coupling%dqsfc_cpl(i) + Coupling%dqsfci_cpl(i) * dtf ! + else + Coupling%dusfc_cpl(i) = huge + Coupling%dvsfc_cpl(i) = huge + Coupling%dtsfc_cpl(i) = huge + Coupling%dqsfc_cpl(i) = huge endif ! Ocean only, NO LAKES enddo endif @@ -3276,10 +3293,10 @@ subroutine GFS_physics_driver & ! print *,' dtdt=',dtdt(ipr,:) ! print *,' gu0=',gu0(ipr,:) ! print *,' gv0=',gv0(ipr,:) -! write(0,*) ' gt0=',(gt0(ipr,k),k=1,levs),' kdt=',kdt -! write(0,*)' gq0=',(gq0(ipr,k,1),k=1,levs),' lat=',lat -! write(0,*)' gq0i2=',(gq0(ipr,k,ntiw),k=1,levs),' lat=',lat -! write(0,*)' gq1=',(gq0(ipr,k,ntcw),k=1,levs) +! write(0,*) ' gt0=',(Stateout%gt0(ipr,k),k=1,levs),' kdt=',kdt +! write(0,*)' gq0=',(Stateout%gq0(ipr,k,1),k=1,levs) +! write(0,*)' gq0i2=',(Stateout%gq0(ipr,k,ntiw),k=1,levs) +! write(0,*)' gq1=',(Stateout%gq0(ipr,k,ntcw),k=1,levs) ! print *,' vvel=',vvel ! endif ! if (lprnt) write(7000,*)' bef convection gu0=',gu0(ipr,:) @@ -3356,6 +3373,7 @@ subroutine GFS_physics_driver & do n=2,ntrac if ( n /= ntcw .and. n /= ntiw .and. n /= ntclamt .and. & n /= ntrw .and. n /= ntsw .and. n /= ntrnc .and. & +! n /= ntlnc .and. n /= ntinc .and. & n /= ntsnc .and. n /= ntgl .and. n /= ntgnc) then tracers = tracers + 1 do k=1,levs @@ -3491,6 +3509,11 @@ subroutine GFS_physics_driver & rhc(:,:) = one !*## CCPP ## endif + +! if (lprnt) write(0,*)' clwice=',clw(ipr,:,1) +! if (lprnt) write(0,*)' clwwat=',clw(ipr,:,2) +! if (lprnt) write(0,*)' rhc=',rhc(ipr,:) + ! ! Call SHOC if do_shoc is true and shocaftcnv is false ! @@ -3556,6 +3579,8 @@ subroutine GFS_physics_driver & ! if (lprnt) write(0,*)'gq01=',Stateout%gq0(ipr,:,1) ! if (lprnt) write(0,*)'clwi=',clw(ipr,:,1) ! if (lprnt) write(0,*)'clwl=',clw(ipr,:,2) +! if (lprnt) write(0,*)'befncpi=',ncpi(ipr,:) +! if (lprnt) write(0,*)'tkh=',Tbd%phy_f3d(ipr,:,ntot3d-1) ! if (lprnt) write(0,*) ' befshoc hflx=',hflx(ipr),' evap=',evap(ipr),& ! ' stress=',stress(ipr) ! dtshoc = 60.0 @@ -3594,6 +3619,7 @@ subroutine GFS_physics_driver & lprnt, ipr, imp_physics, ncpl, ncpi) +! if (lprnt) write(0,*)'aftncpi=',ncpi(ipr,:) ! enddo ! if (imp_physics == Model%imp_physics_mg .and. Model%fprcp > 1) then ! do k=1,levs @@ -3604,7 +3630,7 @@ subroutine GFS_physics_driver & ! endif ! if (lprnt) write(0,*)'aftshocgt0=',Stateout%gt0(ipr,:) -! if (lprnt) write(0,*)'aftshocgq0=',Stateout%gq0(ipr,1:60,1) +! if (lprnt) write(0,*)'aftshocgq0=',Stateout%gq0(ipr,:,1) ! if (lprnt) write(0,*)' aft shoc tke=',clw(ipr,1:25,ntk), & ! &' kdt=',kdt,'xlon=',grid%xlon(ipr),' xlat=',grid%xlat(ipr) ! if (lprnt) write(0,*)' aftshoccld=',tbd%phy_f3d(ipr,:,ntot3d-2)*100 @@ -3918,8 +3944,8 @@ subroutine GFS_physics_driver & trcmin, ntk) !*## CCPP ## -! if (lprnt) write(0,*)' gt04=',Stateout%gt0(ipr,1:60) -! if (lprnt) write(0,*)' gq04=',Stateout%gq0(ipr,1:60,1) +! if (lprnt) write(0,*)' gt04=',Stateout%gt0(ipr,:) +! if (lprnt) write(0,*)' gq04=',Stateout%gq0(ipr,:,1) ! if (lprnt) write(0,*)'aftrasclw1=',clw(ipr,:,1) ! if (lprnt) write(0,*)'aftrasclw2=',clw(ipr,:,2) ! if (lprnt) write(0,*)'aftrastke=',clw(ipr,:,ntk) @@ -4039,6 +4065,13 @@ subroutine GFS_physics_driver & ! !----------------Convective gravity wave drag parameterization starting -------- +! if (lprnt) then +! write(0,*) ' befgwgt0=',Stateout%gt0(ipr,:) +! write(0,*) ' befgwgq0=',Stateout%gq0(ipr,:,1) +! write(0,*) ' do_cnvgwd=',Model%do_cnvgwd +! endif + +! DH* this block is in gwdc_pre !## CCPP ##* gwdc.f/gwdc_pre Note: The conditional above is not in the scheme, so ! the execution of the code below is controlled by its presence in the CCPP SDF ! --- ... calculate maximum convective heating rate @@ -4231,6 +4264,11 @@ subroutine GFS_physics_driver & deallocate(gwdcu, gwdcv) endif ! end if_cnvgwd (convective gravity wave drag) +! if (lprnt) then +! write(0,*) ' befgwegt0=',Stateout%gt0(ipr,:) +! write(0,*) ' befgwegq0=',Stateout%gq0(ipr,:,1) +! endif + ! if (lprnt) write(7000,*)' aft cnvgwd gu0=',gu0(ipr,:) ! if (lprnt) write(7000,*)' aft cnvgwd gv0=',gv0(ipr,:) ! &,' lat=',lat,' kdt=',kdt,' me=',me @@ -4290,6 +4328,10 @@ subroutine GFS_physics_driver & else nsamftrac = tottracer endif +! if (lprnt) then +! write(0,*) ' befshgt0=',Stateout%gt0(ipr,:) +! write(0,*) ' befshgq0=',Stateout%gq0(ipr,:,1) +! endif !*## CCPP ## !## CCPP ##* samfshalcnv.f/samfshalcnv_run call samfshalcnv (im, ix, levs, dtp, itc, Model%ntchm, ntk, nsamftrac, & @@ -4507,6 +4549,7 @@ subroutine GFS_physics_driver & ! if ( n /= ntcw .and. n /= ntiw .and. n /= ntclamt) then if ( n /= ntcw .and. n /= ntiw .and. n /= ntclamt .and. & n /= ntrw .and. n /= ntsw .and. n /= ntrnc .and. & +! n /= ntlnc .and. n /= ntinc .and. & n /= ntsnc .and. n /= ntgl .and. n /= ntgnc ) then tracers = tracers + 1 do k=1,levs @@ -4569,6 +4612,16 @@ subroutine GFS_physics_driver & endif ! end if_ntcw !*## CCPP ## +! if (lprnt) then +! write(0,*)' aft shallow physics kdt=',kdt +! write(0,*)'qt0s=',Stateout%gt0(ipr,:) +! write(0,*)'qq0s=',Stateout%gq0(ipr,:,1) +! write(0,*)'qq0ws=',Stateout%gq0(ipr,:,ntcw) +! write(0,*)'qq0is=',Stateout%gq0(ipr,:,ntiw) +! write(0,*)'qq0ntic=',Stateout%gq0(ipr,:,ntinc) +! write(0,*)'qq0os=',Stateout%gq0(ipr,:,ntoz) +! endif + ! Legacy routine which determines convectve clouds - should be removed at some point !## CCPP ## cnvc90.f/cnvc90_run call cnvc90 (Model%clstp, im, ix, Diag%rainc, kbot, ktop, levs, Statein%prsi, & @@ -4776,12 +4829,13 @@ subroutine GFS_physics_driver & Tbd%phy_f3d(:,:,1),Tbd%phy_f3d(:,:,2),Tbd%phy_f3d(:,:,3), & ims,ime, kms,kme, & its,ite, kts,kte) +! !*## CCPP ## - elseif (imp_physics == Model%imp_physics_mg) then ! MGB double-moment microphysics - ! ------------------------------ + elseif (imp_physics == Model%imp_physics_mg) then ! MGB double-moment microphysics + ! ------------------------------ !## CCPP ##* GFS_typedefs.F90/control_initialize - kk = 5 - if (Model%fprcp >= 2) kk = 6 + kk = 5 + if (Model%fprcp >= 2) kk = 6 !*## CCPP ## ! Acheng used clw here for other code to run smoothly and minimum change ! to make the code work. However, the nc and clw should be treated @@ -4789,86 +4843,86 @@ subroutine GFS_physics_driver & ! year. I believe this will make the physical interaction more reasonable ! Anning 12/5/2015 changed ntcw hold liquid only !## CCPP ##* m_micro_insterstitial.F90/m_micro_pre_run - if (Model%do_shoc) then - skip_macro = Model%do_shoc - if (Model%fprcp == 0) then - do k=1,levs - do i=1,im - clw(i,k,1) = Stateout%gq0(i,k,ntiw) ! ice - clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water - Tbd%phy_f3d(i,k,1) = Tbd%phy_f3d(i,k,ntot3d-2) ! clouds from shoc + if (Model%do_shoc) then + skip_macro = Model%do_shoc + if (Model%fprcp == 0) then + do k=1,levs + do i=1,im + clw(i,k,1) = Stateout%gq0(i,k,ntiw) ! ice + clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water + Tbd%phy_f3d(i,k,1) = Tbd%phy_f3d(i,k,ntot3d-2) ! clouds from shoc + enddo enddo - enddo - elseif (abs(Model%fprcp) == 1 .or. mg3_as_mg2) then - do k=1,levs - do i=1,im - clw(i,k,1) = Stateout%gq0(i,k,ntiw) ! ice - clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water - qrn(i,k) = Stateout%gq0(i,k,ntrw) - qsnw(i,k) = Stateout%gq0(i,k,ntsw) - ncpr(i,k) = Stateout%gq0(i,k,ntrnc) - ncps(i,k) = Stateout%gq0(i,k,ntsnc) - Tbd%phy_f3d(i,k,1) = Tbd%phy_f3d(i,k,ntot3d-2) ! clouds from shoc + elseif (abs(Model%fprcp) == 1 .or. mg3_as_mg2) then + do k=1,levs + do i=1,im + clw(i,k,1) = Stateout%gq0(i,k,ntiw) ! ice + clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water + qrn(i,k) = Stateout%gq0(i,k,ntrw) + qsnw(i,k) = Stateout%gq0(i,k,ntsw) + ncpr(i,k) = Stateout%gq0(i,k,ntrnc) + ncps(i,k) = Stateout%gq0(i,k,ntsnc) + Tbd%phy_f3d(i,k,1) = Tbd%phy_f3d(i,k,ntot3d-2) ! clouds from shoc + enddo enddo - enddo - else - do k=1,levs - do i=1,im - clw(i,k,1) = Stateout%gq0(i,k,ntiw) ! ice - clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water - qrn(i,k) = Stateout%gq0(i,k,ntrw) - qsnw(i,k) = Stateout%gq0(i,k,ntsw) - qgl(i,k) = Stateout%gq0(i,k,ntgl) - ncpr(i,k) = Stateout%gq0(i,k,ntrnc) - ncps(i,k) = Stateout%gq0(i,k,ntsnc) - ncgl(i,k) = Stateout%gq0(i,k,ntgnc) - Tbd%phy_f3d(i,k,1) = Tbd%phy_f3d(i,k,ntot3d-2) ! clouds from shoc + else + do k=1,levs + do i=1,im + clw(i,k,1) = Stateout%gq0(i,k,ntiw) ! ice + clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water + qrn(i,k) = Stateout%gq0(i,k,ntrw) + qsnw(i,k) = Stateout%gq0(i,k,ntsw) + qgl(i,k) = Stateout%gq0(i,k,ntgl) + ncpr(i,k) = Stateout%gq0(i,k,ntrnc) + ncps(i,k) = Stateout%gq0(i,k,ntsnc) + ncgl(i,k) = Stateout%gq0(i,k,ntgnc) + Tbd%phy_f3d(i,k,1) = Tbd%phy_f3d(i,k,ntot3d-2) ! clouds from shoc + enddo enddo - enddo - endif + endif - else + else ! clouds from t-dt and cnvc - if (Model%fprcp == 0 ) then - do k=1,levs - do i=1,im - clw(i,k,1) = Stateout%gq0(i,k,ntiw) ! ice - clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water + if (Model%fprcp == 0 ) then + do k=1,levs + do i=1,im + clw(i,k,1) = Stateout%gq0(i,k,ntiw) ! ice + clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water + enddo enddo - enddo - elseif (abs(Model%fprcp) == 1 .or. mg3_as_mg2) then - do k=1,levs - do i=1,im - clw(i,k,1) = Stateout%gq0(i,k,ntiw) ! ice - clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water - qrn(i,k) = Stateout%gq0(i,k,ntrw) - qsnw(i,k) = Stateout%gq0(i,k,ntsw) - ncpr(i,k) = Stateout%gq0(i,k,ntrnc) - ncps(i,k) = Stateout%gq0(i,k,ntsnc) + elseif (abs(Model%fprcp) == 1 .or. mg3_as_mg2) then + do k=1,levs + do i=1,im + clw(i,k,1) = Stateout%gq0(i,k,ntiw) ! ice + clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water + qrn(i,k) = Stateout%gq0(i,k,ntrw) + qsnw(i,k) = Stateout%gq0(i,k,ntsw) + ncpr(i,k) = Stateout%gq0(i,k,ntrnc) + ncps(i,k) = Stateout%gq0(i,k,ntsnc) + enddo enddo - enddo - else - do k=1,levs - do i=1,im - clw(i,k,1) = Stateout%gq0(i,k,ntiw) ! ice - clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water - qrn(i,k) = Stateout%gq0(i,k,ntrw) - qsnw(i,k) = Stateout%gq0(i,k,ntsw) - qgl(i,k) = Stateout%gq0(i,k,ntgl) - ncpr(i,k) = Stateout%gq0(i,k,ntrnc) - ncps(i,k) = Stateout%gq0(i,k,ntsnc) - ncgl(i,k) = Stateout%gq0(i,k,ntgnc) + else + do k=1,levs + do i=1,im + clw(i,k,1) = Stateout%gq0(i,k,ntiw) ! ice + clw(i,k,2) = Stateout%gq0(i,k,ntcw) ! water + qrn(i,k) = Stateout%gq0(i,k,ntrw) + qsnw(i,k) = Stateout%gq0(i,k,ntsw) + qgl(i,k) = Stateout%gq0(i,k,ntgl) + ncpr(i,k) = Stateout%gq0(i,k,ntrnc) + ncps(i,k) = Stateout%gq0(i,k,ntsnc) + ncgl(i,k) = Stateout%gq0(i,k,ntgnc) + enddo enddo - enddo + endif endif - endif ! add convective cloud fraction - do k = 1,levs - do i = 1,im - Tbd%phy_f3d(i,k,1) = min(one, Tbd%phy_f3d(i,k,1) + clcn(i,k)) + do k = 1,levs + do i = 1,im + Tbd%phy_f3d(i,k,1) = min(one, Tbd%phy_f3d(i,k,1) + clcn(i,k)) + enddo enddo - enddo !*## CCPP ## ! notice clw ix instead of im @@ -4878,6 +4932,7 @@ subroutine GFS_physics_driver & ! if(lprnt) write(0,*) ' befgt0=',Stateout%gt0(ipr,:),' kdt=',kdt ! if(lprnt) write(0,*) ' befgq0=',Stateout%gq0(ipr,:,1),' kdt=',kdt ! if(lprnt) write(0,*) ' befntlnc=',Stateout%gq0(ipr,:,ntlnc),' kdt=',kdt +! if(lprnt) write(0,*) ' befntinc=',Stateout%gq0(ipr,:,ntinc),' kdt=',kdt ! if (lprnt) write(0,*)' clw1bef=',clw(ipr,:,1),' kdt=',kdt ! if (lprnt) write(0,*)' clw2bef=',clw(ipr,:,2),' kdt=',kdt ! if (lprnt) write(0,*)' qrnb=',qrn(ipr,:),' kdt=',kdt @@ -4891,31 +4946,32 @@ subroutine GFS_physics_driver & ! do k=1,levs ! write(1000+me,*)' maxwatncb=',maxval(Stateout%gq0(1:im,k,ntlnc)),' k=',k,' kdt',kdt ! enddo + !## CCPP ##* m_micro.F90/m_micro_run - call m_micro_driver (im, ix, levs, Model%flipv, dtp, Statein%prsl, & - Statein%prsi, Statein%phil, Statein%phii, & - Statein%vvl, clw(1,1,2), QLCN, clw(1,1,1), QICN, & - Radtend%htrlw, Radtend%htrsw, w_upi, cf_upi, & - FRLAND, Diag%HPBL, CNV_MFD, CNV_DQLDT, & -! FRLAND, Diag%HPBL, CNV_MFD, CNV_PRC3, CNV_DQLDT, & - CLCN, Stateout%gu0, Stateout%gv0, Diag%dusfc, & - Diag%dvsfc, dusfc1, dvsfc1, dusfc1, dvsfc1, & - CNV_FICE, CNV_NDROP, CNV_NICE, Stateout%gq0(1,1,1), & - Stateout%gq0(1,1,ntcw), & - Stateout%gq0(1,1,ntiw), Stateout%gt0, rain1, & - Diag%sr, Stateout%gq0(1,1,ntlnc), & - Stateout%gq0(1,1,ntinc), Model%fprcp, qrn, & - qsnw, qgl, ncpr, ncps, ncgl, & - Tbd%phy_f3d(1,1,1), kbot, & - Tbd%phy_f3d(1,1,2), Tbd%phy_f3d(1,1,3), & - Tbd%phy_f3d(1,1,4), Tbd%phy_f3d(1,1,5), & - Tbd%phy_f3d(1,1,kk), Tbd%aer_nm, & - Model%aero_in, Tbd%in_nm, Tbd%ccn_nm, Model%iccn, & - skip_macro, lprnt, & -! skip_macro, cn_prc, cn_snr, lprnt, & -! ipr, kdt, Grid%xlat, Grid%xlon) - Model%mg_alf, Model%mg_qcmin, Model%pdfflag, & - ipr, kdt, Grid%xlat, Grid%xlon, rhc) + call m_micro_driver (im, ix, levs, Model%flipv, dtp, Statein%prsl, & + Statein%prsi, Statein%phil, Statein%phii, & + Statein%vvl, clw(1,1,2), QLCN, clw(1,1,1), QICN, & + Radtend%htrlw, Radtend%htrsw, w_upi, cf_upi, & + FRLAND, Diag%HPBL, CNV_MFD, CNV_DQLDT, & +! FRLAND, Diag%HPBL, CNV_MFD, CNV_PRC3, CNV_DQLDT, & + CLCN, Stateout%gu0, Stateout%gv0, Diag%dusfc, & + Diag%dvsfc, dusfc1, dvsfc1, dusfc1, dvsfc1, & + CNV_FICE, CNV_NDROP, CNV_NICE, Stateout%gq0(1,1,1), & + Stateout%gq0(1,1,ntcw), & + Stateout%gq0(1,1,ntiw), Stateout%gt0, rain1, & + Diag%sr, Stateout%gq0(1,1,ntlnc), & + Stateout%gq0(1,1,ntinc), Model%fprcp, qrn, & + qsnw, qgl, ncpr, ncps, ncgl, & + Tbd%phy_f3d(1,1,1), kbot, & + Tbd%phy_f3d(1,1,2), Tbd%phy_f3d(1,1,3), & + Tbd%phy_f3d(1,1,4), Tbd%phy_f3d(1,1,5), & + Tbd%phy_f3d(1,1,kk), Tbd%aer_nm, & + Model%aero_in, Tbd%in_nm, Tbd%ccn_nm, Model%iccn, & + skip_macro, lprnt, & +! skip_macro, cn_prc, cn_snr, lprnt, & +! ipr, kdt, Grid%xlat, Grid%xlon) + Model%mg_alf, Model%mg_qcmin, Model%pdfflag, & + ipr, kdt, Grid%xlat, Grid%xlon, rhc) !*## CCPP ## ! do k=1,levs ! write(1000+me,*)' maxwatnca=',maxval(Stateout%gq0(1:im,k,ntlnc)),' k=',k,' kdt=',kdt @@ -4936,7 +4992,7 @@ subroutine GFS_physics_driver & ! &,' cn_prc=',cn_prc(ipr),' cn_snr=',cn_snr(ipr),' kdt=',kdt ! if(lprnt) write(0,*) ' aftgt0=',Stateout%gt0(ipr,:),' kdt=',kdt ! if (lprnt) write(0,*) ' aftlsgq0=',stateout%gq0(ipr,:,1),' kdt=',kdt -! if (lprnt) write(0,*)' clw1aft=',stateout%gq0(ipr,:,ntiw),' kdt=',kdt +! if (lprnt) write(0,*)' cli1aft=',stateout%gq0(ipr,:,ntiw),' kdt=',kdt ! if (ntgl > 0 .and. lprnt) & ! write(0,*)' cgw1aft=',stateout%gq0(ipr,:,ntgl),' kdt=',kdt ! if (lprnt) write(0,*)' cloudsm=',tbd%phy_f3d(ipr,:,1)*100,' kdt=',kdt @@ -4944,44 +5000,47 @@ subroutine GFS_physics_driver & ! if (lprnt) write(0,*)' qrna=',qrn(ipr,:),' kdt=',kdt ! if (lprnt) write(0,*)' qsnwa=',qsnw(ipr,:),' kdt=',kdt ! if (lprnt) write(0,*)' qglba',qgl(ipr,:),' kdt=',kdt + !## CCPP ##* m_micro_interstitial.F90/m_micro_post_run - tem = dtp * con_p001 / con_day - if (abs(Model%fprcp) == 1 .or. mg3_as_mg2) then - do k=1,levs + + tem = dtp * con_p001 / con_day + if (abs(Model%fprcp) == 1 .or. mg3_as_mg2) then + do k=1,levs + do i=1,im + if (abs(qrn(i,k)) < qsmall) qrn(i,k) = zero + if (abs(qsnw(i,k)) < qsmall) qsnw(i,k) = zero + Stateout%gq0(i,k,ntrw) = qrn(i,k) + Stateout%gq0(i,k,ntsw) = qsnw(i,k) + Stateout%gq0(i,k,ntrnc) = ncpr(i,k) + Stateout%gq0(i,k,ntsnc) = ncps(i,k) + enddo + enddo do i=1,im - if (abs(qrn(i,k)) < qsmall) qrn(i,k) = zero - if (abs(qsnw(i,k)) < qsmall) qsnw(i,k) = zero - Stateout%gq0(i,k,ntrw) = qrn(i,k) - Stateout%gq0(i,k,ntsw) = qsnw(i,k) - Stateout%gq0(i,k,ntrnc) = ncpr(i,k) - Stateout%gq0(i,k,ntsnc) = ncps(i,k) + Diag%ice(i) = tem * Stateout%gq0(i,1,ntiw) + Diag%snow(i) = tem * qsnw(i,1) + enddo + elseif (Model%fprcp > 1) then + do k=1,levs + do i=1,im + if (abs(qrn(i,k)) < qsmall) qrn(i,k) = zero + if (abs(qsnw(i,k)) < qsmall) qsnw(i,k) = zero + if (abs(qgl(i,k)) < qsmall) qgl(i,k) = zero + Stateout%gq0(i,k,ntrw) = qrn(i,k) + Stateout%gq0(i,k,ntsw) = qsnw(i,k) + Stateout%gq0(i,k,ntgl) = qgl(i,k) + Stateout%gq0(i,k,ntrnc) = ncpr(i,k) + Stateout%gq0(i,k,ntsnc) = ncps(i,k) + Stateout%gq0(i,k,ntgnc) = ncgl(i,k) + enddo enddo - enddo - do i=1,im - Diag%ice(i) = tem * Stateout%gq0(i,1,ntiw) - Diag%snow(i) = tem * qsnw(i,1) - enddo - elseif (Model%fprcp > 1) then - do k=1,levs do i=1,im - if (abs(qrn(i,k)) < qsmall) qrn(i,k) = zero - if (abs(qsnw(i,k)) < qsmall) qsnw(i,k) = zero - if (abs(qgl(i,k)) < qsmall) qgl(i,k) = zero - Stateout%gq0(i,k,ntrw) = qrn(i,k) - Stateout%gq0(i,k,ntsw) = qsnw(i,k) - Stateout%gq0(i,k,ntgl) = qgl(i,k) - Stateout%gq0(i,k,ntrnc) = ncpr(i,k) - Stateout%gq0(i,k,ntsnc) = ncps(i,k) - Stateout%gq0(i,k,ntgnc) = ncgl(i,k) + Diag%ice(i) = tem * Stateout%gq0(i,1,ntiw) + Diag%snow(i) = tem * qsnw(i,1) + Diag%graupel(i) = tem * qgl(i,1) enddo - enddo - do i=1,im - Diag%ice(i) = tem * Stateout%gq0(i,1,ntiw) - Diag%snow(i) = tem * qsnw(i,1) - Diag%graupel(i) = tem * qgl(i,1) - enddo + + endif !*## CCPP ## - endif ! if (lprnt) write(0,*)' cloudsm=',tbd%phy_f3d(ipr,:,1)*100,' kdt=',kdt ! if (lprnt) write(0,*)' clw2aft=',stateout%gq0(ipr,:,ntcw),' kdt=',kdt @@ -5050,20 +5109,20 @@ subroutine GFS_physics_driver & reset) tem = dtp * con_p001 / con_day do i = 1, im -! rain0(i,1) = max(con_d00, rain0(i,1)) -! snow0(i,1) = max(con_d00, snow0(i,1)) -! ice0(i,1) = max(con_d00, ice0(i,1)) -! graupel0(i,1) = max(con_d00, graupel0(i,1)) - if(rain0(i,1)*tem < rainmin) then - rain0(i,1) = zero +! rain0(i,1) = max(con_d00, rain0(i,1)) +! snow0(i,1) = max(con_d00, snow0(i,1)) +! ice0(i,1) = max(con_d00, ice0(i,1)) +! graupel0(i,1) = max(con_d00, graupel0(i,1)) + if (rain0(i,1)*tem < rainmin) then + rain0(i,1) = zero endif - if(ice0(i,1)*tem < rainmin) then + if (ice0(i,1)*tem < rainmin) then ice0(i,1) = zero endif - if(snow0(i,1)*tem < rainmin) then + if (snow0(i,1)*tem < rainmin) then snow0(i,1) = zero endif - if(graupel0(i,1)*tem < rainmin) then + if (graupel0(i,1)*tem < rainmin) then graupel0(i,1) = zero endif @@ -5103,7 +5162,7 @@ subroutine GFS_physics_driver & enddo - if(Model%effr_in) then + if (Model%effr_in) then do i =1, im den(i,k) = 0.622*Statein%prsl(i,k) / & (con_rd*Stateout%gt0(i,k)*(Stateout%gq0(i,k,1)+0.622)) @@ -5113,25 +5172,25 @@ subroutine GFS_physics_driver & !*## CCPP ## !## CCPP ##* maximum_hourly_diagnostics.F90/maximum_hourly_diagnsostics_run !Calculate hourly max 1-km agl and -10C reflectivity - if (Model%lradar .and. & - (imp_physics == Model%imp_physics_gfdl .or. & - imp_physics == Model%imp_physics_thompson)) then - allocate(refd(im)) - allocate(refd263k(im)) - call max_fields(Statein%phil,Diag%refl_10cm,con_g,im,levs,refd,Stateout%gt0,refd263k) - if (reset) then + if (Model%lradar .and. & + (imp_physics == Model%imp_physics_gfdl .or. & + imp_physics == Model%imp_physics_thompson)) then + allocate(refd(im)) + allocate(refd263k(im)) + call max_fields(Statein%phil,Diag%refl_10cm,con_g,im,levs,refd,Stateout%gt0,refd263k) + if (reset) then + do i=1,im + Diag%refdmax(I) = -35. + Diag%refdmax263k(I) = -35. + enddo + endif do i=1,im - Diag%refdmax(I) = -35. - Diag%refdmax263k(I) = -35. + Diag%refdmax(i) = max(Diag%refdmax(i),refd(i)) + Diag%refdmax263k(i) = max(Diag%refdmax263k(i),refd263k(i)) enddo + deallocate (refd) + deallocate (refd263k) endif - do i=1,im - Diag%refdmax(i) = max(Diag%refdmax(i),refd(i)) - Diag%refdmax263k(i) = max(Diag%refdmax263k(i),refd263k(i)) - enddo - deallocate (refd) - deallocate (refd263k) - endif !*## CCPP ## !## CCPP ##* gfdl_cloud_microphys.F90/gfdl_cloud_microphys_run if(Model%effr_in) then @@ -5144,23 +5203,23 @@ subroutine GFS_physics_driver & Tbd%phy_f3d(1:im,1:levs,1), Tbd%phy_f3d(1:im,1:levs,2), & Tbd%phy_f3d(1:im,1:levs,3), Tbd%phy_f3d(1:im,1:levs,4), & Tbd%phy_f3d(1:im,1:levs,5)) + !*## CCPP ## -! do k = 1, levs -! do i=1,im -! -! if(Model%me==0) then -! if(Tbd%phy_f3d(i,k,1) > 5.) then -! write(6,*) 'phy driver:cloud radii:',Model%kdt, i,k, & -! Tbd%phy_f3d(i,k,1) -! endif -! if(Tbd%phy_f3d(i,k,3)> zero) then -! write(6,*) 'phy driver:rain radii:',Model%kdt, i,k, & -! Tbd%phy_f3d(i,k,3) -! endif -! -! endif -! enddo -! enddo +! do k = 1, levs +! do i=1,im +! if(Model%me==0) then +! if(Tbd%phy_f3d(i,k,1) > 5.) then +! write(6,*) 'phy driver:cloud radii:',Model%kdt, i,k, & +! Tbd%phy_f3d(i,k,1) +! endif +! if(Tbd%phy_f3d(i,k,3)> zero) then +! write(6,*) 'phy driver:rain radii:',Model%kdt, i,k, & +! Tbd%phy_f3d(i,k,3) +! endif +! +! endif +! enddo +! enddo endif @@ -5220,9 +5279,10 @@ subroutine GFS_physics_driver & rain1(i) = max(rain1(i) - temrain1(i)*0.001, 0.0_kind_phys) enddo endif + !*## CCPP ## !## CCPP ##* GFS_MP_generic.F90/GFS_MP_generic_post_run - Diag%rain(:) = Diag%rainc(:) + frain * rain1(:) + Diag%rain(:) = Diag%rainc(:) + frain * rain1(:) ! total rain per timestep ! --- get the amount of different precip type for Noah MP ! --- convert from m/dtp to mm/s @@ -5378,6 +5438,7 @@ subroutine GFS_physics_driver & enddo elseif( .not. Model%cal_pre) then if (Model%imp_physics == Model%imp_physics_mg) then ! MG microphysics + tem = con_day / (dtp * con_p001) ! mm / day do i=1,im Sfcprop%tprcp(i) = max(zero, Diag%rain(i) ) ! clu: rain -> tprcp if (Diag%rain(i)*tem > rainmin) then @@ -5398,7 +5459,6 @@ subroutine GFS_physics_driver & endif - ! --- ... coupling insertion if (Model%cplflx .or. Model%cplchm) then @@ -5536,7 +5596,15 @@ subroutine GFS_physics_driver & ! &' rain=',rain(ipr),' rainc=',rainc(ipr) ! if (lprnt) call mpi_quit(7) ! if (kdt > 2 ) call mpi_quit(70) -! if (lprnt) write(0,*)'qt0out=',Stateout%gt0(ipr,:) & +! if (lprnt) then +! write(0,*)' at the end of physics kdt=',kdt +! write(0,*)' end rain=',diag%rain(ipr),' rainc=',diag%rainc(ipr) +! write(0,*)'qt0out=',Stateout%gt0(ipr,:) +! write(0,*)'qq0outv=',Stateout%gq0(ipr,:,1) +! write(0,*)'qq0outw=',Stateout%gq0(ipr,:,ntcw) +! write(0,*)'qq0outi=',Stateout%gq0(ipr,:,ntiw) +! write(0,*)'qq0outo=',Stateout%gq0(ipr,:,ntoz) +! endif ! if (lprnt) write(0,*)'gq0outtke=',Stateout%gq0(ipr,1:25,ntke) & ! ,'xlon=',grid%xlon(ipr)*rad2dg,' xlat=',grid%xlat(ipr)*rad2dg ! if (lprnt) write(0,*)' clouddriverend=',Tbd%phy_f3d(ipr,:,1)*100,' kdt=',kdt diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index cf2589755..652955816 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -627,6 +627,8 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) if (nint(oro_var2(1,1,18)) == -9999._kind_phys) then ! lakefrac doesn't exist in the restart, need to create it if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - will computing lakefrac') Model%frac_grid = .false. + elseif (Model%frac_grid_off) then + Model%frac_grid = .false. else Model%frac_grid = .true. endif @@ -936,7 +938,13 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) Sfcprop(nb)%zorll(ix) = sfc_var2(i,j,34) !--- zorll (zorl on land portion of a cell) end if - if(Model%frac_grid) then ! obtain slmsk from landfrac and fice + if(Model%frac_grid) then ! obtain slmsk from landfrac +!! next 5 lines are temporary till lake model is available + if (Sfcprop(nb)%lakefrac(ix) > 0.0) then + Sfcprop(nb)%lakefrac(ix) = nint(Sfcprop(nb)%lakefrac(ix)) + Sfcprop(nb)%landfrac(ix) = 1.-Sfcprop(nb)%lakefrac(ix) + if (Sfcprop(nb)%lakefrac(ix) == 0) Sfcprop(nb)%fice(ix)=0. + end if Sfcprop(nb)%slmsk(ix) = ceiling(Sfcprop(nb)%landfrac(ix)) if (Sfcprop(nb)%fice(ix) > 0. .and. Sfcprop(nb)%landfrac(ix)==0.) Sfcprop(nb)%slmsk(ix) = 2 ! land dominates ice if co-exist else ! obtain landfrac from slmsk From b8f6e8e3a816956771077dae39408e3cf10502e9 Mon Sep 17 00:00:00 2001 From: "Shan.Sun" Date: Thu, 12 Mar 2020 19:18:55 +0000 Subject: [PATCH 12/15] Setting fice=1 if fice>1.-1.e-10 --- atmos_model.F90 | 1 + gfsphysics/GFS_layer/GFS_physics_driver.F90 | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/atmos_model.F90 b/atmos_model.F90 index b9a91b94c..607112e19 100644 --- a/atmos_model.F90 +++ b/atmos_model.F90 @@ -1747,6 +1747,7 @@ subroutine assign_importdata(rc) IPD_Data(nb)%Coupling%slimskin_cpl(ix) = IPD_Data(nb)%Sfcprop%slmsk(ix) if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then IPD_Data(nb)%Coupling%ficein_cpl(ix) = max(zero, min(one, datar8(i,j)/IPD_Data(nb)%Sfcprop%oceanfrac(ix))) !LHS: ice frac wrt water area + if (IPD_Data(nb)%Coupling%ficein_cpl(ix) > one-epsln) IPD_Data(nb)%Coupling%ficein_cpl(ix)=one if (IPD_Data(nb)%Coupling%ficein_cpl(ix) >= IPD_control%min_seaice) then if (abs(one-IPD_Data(nb)%Sfcprop%oceanfrac(ix)) < epsln) IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points IPD_Data(nb)%Coupling%slimskin_cpl(ix) = 4. diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index 22bc0e832..b1f0c7d1c 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -2838,7 +2838,7 @@ subroutine GFS_physics_driver & if (Model%cplflx) then do i=1,im if (Sfcprop%oceanfrac(i) > zero) then ! Ocean only, NO LAKES - if (Sfcprop%fice(i) >= one - epsln) then ! no open water, thus use results from CICE + if (Sfcprop%fice(i) > one - epsln) then ! no open water, thus use results from CICE Coupling%dusfci_cpl(i) = dusfc_cice(i) Coupling%dvsfci_cpl(i) = dvsfc_cice(i) Coupling%dtsfci_cpl(i) = dtsfc_cice(i) From 8645684f2d3b341796c435977e247f339ea8d873 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Sun, 15 Mar 2020 19:30:08 -0600 Subject: [PATCH 13/15] gfsphysics/GFS_layer/GFS_typedefs.F90: bugfix for Thompson radar reflectivity reset logic --- gfsphysics/GFS_layer/GFS_typedefs.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gfsphysics/GFS_layer/GFS_typedefs.F90 b/gfsphysics/GFS_layer/GFS_typedefs.F90 index ba3cc58f7..8250bea53 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -6528,7 +6528,7 @@ subroutine interstitial_phys_reset (Interstitial, Model) ! Set flag for resetting maximum hourly output fields Interstitial%reset = mod(Model%kdt-1, nint(Model%avg_max_length/Model%dtp)) == 0 ! Set flag for resetting radar reflectivity calculation - if (Interstitial%radar_reset<0) then + if (Model%nsradar_reset<0) then Interstitial%radar_reset = .true. else Interstitial%radar_reset = mod(Model%kdt-1, nint(Model%nsradar_reset/Model%dtp)) == 0 From 811ca4da4ad3d3c481f0a8b21331f745e51c04d5 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Sun, 15 Mar 2020 19:30:22 -0600 Subject: [PATCH 14/15] Support for CCPP with hera.gnu --- ccpp/build_ccpp.sh | 8 +++++++- ccpp/set_compilers.sh | 7 +++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/ccpp/build_ccpp.sh b/ccpp/build_ccpp.sh index f6bec6ff9..b5b8c53d7 100755 --- a/ccpp/build_ccpp.sh +++ b/ccpp/build_ccpp.sh @@ -5,7 +5,7 @@ set -eu # List of valid/tested machines VALID_MACHINES=( wcoss_cray wcoss_dell_p3 gaea.intel jet.intel \ - hera.intel \ + hera.intel hera.gnu \ cheyenne.intel cheyenne.intel-impi cheyenne.gnu cheyenne.pgi endeavor.intel \ stampede.intel supermuc_phase2.intel macosx.gnu \ linux.intel linux.gnu linux.pgi ) @@ -128,6 +128,12 @@ fi if [[ "${MAKE_OPT}" == *"STATIC=Y"* ]]; then CCPP_CMAKE_FLAGS="${CCPP_CMAKE_FLAGS} -DSTATIC=ON" else + # hera.gnu uses the NCEPLIBS-external/NCEPLIBS umbrella build libraries, + # which cannot be linked dynamically at this point (missing -fPIC flag) + if [[ "${MACHINE_ID}" == "hera.gnu" ]]; then + echo "Dynamic CCPP build not supported on hera.gnu at this time." + exit 1 + fi # Dynamic builds require linking the NCEPlibs, provide path to them CCPP_CMAKE_FLAGS="${CCPP_CMAKE_FLAGS} -DSTATIC=OFF -DBACIO_LIB4=${BACIO_LIB4} -DSP_LIBd=${SP_LIBd} -DW3NCO_LIBd=${W3NCO_LIBd}" fi diff --git a/ccpp/set_compilers.sh b/ccpp/set_compilers.sh index 0f591d6c0..5d797a161 100755 --- a/ccpp/set_compilers.sh +++ b/ccpp/set_compilers.sh @@ -39,6 +39,13 @@ case "$MACHINE_ID" in export F77=mpiifort export F90=mpiifort ;; + hera.gnu) + export CC=mpicc + export CXX=mpicxx + export FC=mpif90 + export F77=mpif77 + export F90=mpif90 + ;; cheyenne.intel) export CC=mpicc export CXX=mpicxx From 9b79ce1797cd9564ab097f3e3cc8042562c36ac6 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 17 Mar 2020 08:07:28 -0600 Subject: [PATCH 15/15] Revert change to .gitmodules and update submodule pointer for ccpp-physics --- .gitmodules | 6 ++---- ccpp/physics | 2 +- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/.gitmodules b/.gitmodules index 89bf9d1be..2fdeca40a 100644 --- a/.gitmodules +++ b/.gitmodules @@ -8,7 +8,5 @@ branch = dtc/develop [submodule "ccpp/physics"] path = ccpp/physics - #url = https://github.com/NCAR/ccpp-physics - #branch = dtc/develop - url = https://github.com/climbfuji/ccpp-physics - branch = final_pr_before_merging_to_develop_or_master_20200313 + url = https://github.com/NCAR/ccpp-physics + branch = dtc/develop diff --git a/ccpp/physics b/ccpp/physics index 215399e4c..82012f608 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit 215399e4cfc08855e1d98d5b52cad0967eb93920 +Subproject commit 82012f60829997e18f7bca28e546f8b426e00764