@@ -35,7 +35,7 @@ module satmedmfvdifq
35
35
!!
36
36
!! Incorporate the TTE- EDMF; if (tte_edmf= .true. ),
37
37
!! TKE- EDMF scheme becomes TTE- EDMF scheme and the variable ' te'
38
- !! is read as TTE; if (tte_edmf= .false. ), the variable ' te' is
38
+ !! is read as TTE; if (tte_edmf= .false. ), the variable ' te' is
39
39
!! read as TKE, 5 / 22 / 2025
40
40
!!
41
41
!!
@@ -85,25 +85,25 @@ end subroutine satmedmfvdifq_init
85
85
!! (mfpbltq.f), is represented using a mass flux approach (Siebesma et al.(2007 ) \cite Siebesma_2007 ).
86
86
!! - # A mass- flux approach is also used to represent the stratocumulus- top- induced turbulence
87
87
!! (mfscuq.f).
88
- !! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm
88
+ !! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm
89
89
subroutine satmedmfvdifq_run (im ,km ,ntrac ,ntcw ,ntrw , &
90
90
& ntiw ,ntke ,grav ,pi ,rd ,cp ,rv ,hvap ,hfus ,fv ,eps ,epsm1 , &
91
91
!The following three variables are for SA-3D-TKE
92
92
& def_1 ,def_2 ,def_3 ,sa3dtke ,dku3d_h ,dku3d_e , &
93
- & dv ,du ,tdt ,rtg ,u1 ,v1 ,t1 ,q1 ,usfco ,vsfco ,icplocn2atm , &
93
+ & dv ,du ,tdt ,rtg ,u1 ,v1 ,t1 ,q1 ,usfco ,vsfco ,use_oceanuv , &
94
94
& swh ,hlw ,xmu ,garea ,zvfun ,sigmaf , &
95
95
& psk ,rbsoil ,zorl ,u10m ,v10m ,fm ,fh , &
96
96
& tsea ,heat ,evap ,stress ,spd1 ,kpbl , &
97
97
& prsi ,del ,prsl ,prslk ,phii ,phil ,delt ,tte_edmf , &
98
98
& dspheat ,dusfc ,dvsfc ,dtsfc ,dqsfc ,hpbl ,dkt ,dku ,tkeh , &
99
99
& kinver ,xkzm_m ,xkzm_h ,xkzm_s ,dspfac ,bl_upfr ,bl_dnfr , &
100
100
& rlmx ,elmx ,sfc_rlm ,tc_pbl ,use_lpt , &
101
- !IVAI: canopy inputs from AQM
101
+ !IVAI: canopy inputs from AQM
102
102
& do_canopy , cplaqm , claie , cfch , cfrt , cclu , cpopu , &
103
103
!IVAI
104
104
& ntqv ,dtend ,dtidx ,index_of_temperature ,index_of_x_wind , &
105
105
& index_of_y_wind ,index_of_process_pbl ,gen_tend ,ldiag3d , &
106
- & errmsg ,errflg )
106
+ & errmsg ,errflg )
107
107
!IVAI: aux arrays
108
108
! & naux2d,naux3d,aux2d,aux3d)
109
109
@@ -158,7 +158,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
158
158
& dtend
159
159
integer , intent (in ) :: dtidx(:,:), index_of_temperature, &
160
160
& index_of_x_wind, index_of_y_wind, index_of_process_pbl
161
- integer , intent (in ) :: icplocn2atm
161
+ logical , intent (in ) :: use_oceanuv
162
162
real (kind= kind_phys), intent (out ) :: &
163
163
& dusfc(:), dvsfc(:), &
164
164
& dtsfc(:), dqsfc(:), &
@@ -192,7 +192,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
192
192
integer lcld(im),kcld(im),krad(im),mrad(im)
193
193
integer kx1(im), kb1(im), kpblx(im)
194
194
!
195
- real (kind= kind_phys) te(im,km), tei(im,km-1 ), tke(im,km),
195
+ real (kind= kind_phys) te(im,km), tei(im,km-1 ), tke(im,km),
196
196
& tteh(im,km), tesq(im,km-1 ),e2(im,0 :km)
197
197
!
198
198
real (kind= kind_phys) theta(im,km),thvx(im,km), thlvx(im,km),
@@ -227,7 +227,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
227
227
!
228
228
real (kind= kind_phys) elm(im,km), ele(im,km),
229
229
& ckz(im,km), chz(im,km),
230
- & diss(im,km-1 ),prod(im,km-1 ),
230
+ & diss(im,km-1 ),prod(im,km-1 ),
231
231
& bf(im,km-1 ), shr2(im,km-1 ), wush(im,km),
232
232
& xlamue(im,km-1 ), xlamde(im,km-1 ),
233
233
& gotvx(im,km), rlam(im,km-1 )
@@ -317,7 +317,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
317
317
real (kind= kind_phys) h1
318
318
319
319
real (kind= kind_phys) bfac, mffac
320
-
320
+
321
321
real (kind= kind_phys) qice(im,km),qliq(im,km)
322
322
323
323
!PCC_CANOPY------------------------------------
@@ -438,7 +438,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
438
438
km1 = km - 1
439
439
kmpbl = km / 2
440
440
kmscu = km / 2
441
- !> - Compute physical height of the layer centers and interfaces from
441
+ !> - Compute physical height of the layer centers and interfaces from
442
442
!! the geopotential height (\p zi and \p zl)
443
443
do k= 1 ,km
444
444
do i= 1 ,im
@@ -466,7 +466,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
466
466
do i= 1 ,im
467
467
gdx(i) = sqrt (garea(i))
468
468
enddo
469
- !> - Initialize tke value at vertical layer centers and interfaces
469
+ !> - Initialize tke value at vertical layer centers and interfaces
470
470
!! from tracer (\p tke and \p tkeh)
471
471
do k= 1 ,km
472
472
do i= 1 ,im
@@ -583,11 +583,11 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
583
583
kcld(i) = km1
584
584
endif
585
585
enddo
586
- !
586
+ !
587
587
!> - Compute a function for green vegetation fraction and surface roughness.
588
588
!! Entrainment rate in updraft is a function of vegetation fraction and surface
589
589
!! roughness length
590
- !
590
+ !
591
591
do i = 1 ,im
592
592
tem = (sigmaf(i) - vegflo) / (vegfup - vegflo)
593
593
tem = min (max (tem, 0 .), 1 .)
@@ -598,7 +598,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
598
598
enddo
599
599
!
600
600
!> - Compute \f$\theta\f$(theta), and \f$q_l\f$(qlx), \f$\theta_e\f$(thetae),
601
- !! \f$\theta_v\f$(thvx),\f$\theta_{l,v}\f$ (thlvx) including ice water
601
+ !! \f$\theta_v\f$(thvx),\f$\theta_{l,v}\f$ (thlvx) including ice water
602
602
do k= 1 ,km
603
603
do i= 1 ,im
604
604
pix(i,k) = psk(i) / prslk(i,k)
@@ -637,7 +637,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
637
637
enddo
638
638
enddo
639
639
!
640
- !> - Compute an empirical cloud fraction based on
640
+ !> - Compute an empirical cloud fraction based on
641
641
!! Xu and Randall (1996 ) \cite xu_and_randall_1996
642
642
do k = 1 , km
643
643
do i = 1 , im
@@ -688,7 +688,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
688
688
!
689
689
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
690
690
!
691
- !> - Initialize diffusion coefficients to 0 and calculate the total
691
+ !> - Initialize diffusion coefficients to 0 and calculate the total
692
692
!! radiative heating rate (dku, dkt, radx)
693
693
do k= 1 ,km
694
694
do i= 1 ,im
@@ -705,7 +705,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
705
705
enddo
706
706
!> - Compute stable/ unstable PBL flag (pblflg) based on the total
707
707
!! surface energy flux (\e false if the total surface energy flux
708
- !! is into the surface)
708
+ !! is into the surface)
709
709
do i = 1 ,im
710
710
sflux(i) = heat(i) + evap(i)* fv* theta(i,1 )
711
711
if (.not. sfcflg(i) .or. sflux(i) <= 0 .) pblflg(i)= .false.
@@ -716,12 +716,12 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
716
716
!! - Compute critical bulk Richardson number (\f$Rb_{cr}\f$) (crb)
717
717
!! - For the unstable PBL, crb is a constant (0.25 )
718
718
!! - For the stable boundary layer (SBL), \f$Rb_{cr}\f$ varies
719
- !! with the surface Rossby number, \f$R_{0 }\f$, as given by
719
+ !! with the surface Rossby number, \f$R_{0 }\f$, as given by
720
720
!! Vickers and Mahrt (2004 ) \cite Vickers_2004
721
721
!! \f[
722
722
!! Rb_{cr}= 0.16 (10 ^{- 7 }R_{0 })^{- 0.18 }
723
723
!! \f]
724
- !! \f[
724
+ !! \f[
725
725
!! R_{0 }= \frac{U_{10 }}{f_{0 }z_{0 }}
726
726
!! \f]
727
727
!! where \f$U_{10 }\f$ is the wind speed at 10m above the ground surface,
@@ -753,7 +753,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
753
753
ustar(i) = sqrt (stress(i))
754
754
enddo
755
755
!
756
- !> - Compute buoyancy \f$\frac{\partial \theta_v}{\partial z}\f$ (bf)
756
+ !> - Compute buoyancy \f$\frac{\partial \theta_v}{\partial z}\f$ (bf)
757
757
!! and the wind shear squared (shr2)
758
758
!
759
759
do k = 1 , km1
@@ -980,7 +980,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
980
980
!
981
981
!> - The \f$z/L\f$ (zol) is used as the stability criterion for the PBL.Currently,
982
982
!! strong unstable (convective) PBL for \f$z/L < -0.02\f$ and weakly and moderately
983
- !! unstable PBL for \f$0>z/L>-0.02\f$
983
+ !! unstable PBL for \f$0>z/L>-0.02\f$
984
984
!> - Compute the velocity scale \f$w_s\f$ (wscale) (eqn 22 of Han et al. 2019). It
985
985
!! is represented by the value scaled at the top of the surface layer:
986
986
!! \f[
@@ -1172,11 +1172,11 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
1172
1172
do i = 1, im
1173
1173
if(scuflg(i) .and. kcld(i)==km1) scuflg(i)=.false.
1174
1174
enddo
1175
- !> - Starting at the PBL top and going downward, if the level is less
1175
+ !> - Starting at the PBL top and going downward, if the level is less
1176
1176
!! than the cloud top, find the level of the minimum radiative heating
1177
1177
!! rate wihin the cloud. If the level of the minimum is the lowest model
1178
1178
!! level or the minimum radiative heating rate is positive, then set
1179
- !! scuflg to F.
1179
+ !! scuflg to F.
1180
1180
do i = 1, im
1181
1181
flg(i)=scuflg(i)
1182
1182
enddo
@@ -1202,7 +1202,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
1202
1202
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1203
1203
!> ## Compute components for mass flux mixing by large thermals
1204
1204
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1205
- !> - If the PBL is convective, the updraft properties are initialized
1205
+ !> - If the PBL is convective, the updraft properties are initialized
1206
1206
!! to be the same as the state variables.
1207
1207
do k = 1, km
1208
1208
do i = 1, im
@@ -1236,7 +1236,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
1236
1236
& pcnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx,
1237
1237
& gdx,hpbl,kpbl,vpert,buou,wush,tkemean,vez0fun,xmf,
1238
1238
& tcko,qcko,ucko,vcko,xlamue,bl_upfr)
1239
- !> - Call mfscuq(), which is a new mass-flux parameterization for
1239
+ !> - Call mfscuq(), which is a new mass-flux parameterization for
1240
1240
!! stratocumulus-top-induced turbulence mixing. For details of the mfscuq subroutine, step into its documentation ::mfscuq
1241
1241
call mfscuq(im,im,km,kmscu,ntcw,ntrac1,dt2,
1242
1242
& scuflg,zl,zm,q1,t1,u1,v1,plyr,pix,
@@ -1389,7 +1389,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
1389
1389
!! l_2=min(l_{up},l_{down})
1390
1390
!!\f]
1391
1391
!! and dissipation length scale \f$l_d\f$ is given by:
1392
- !!\f[
1392
+ !!\f[
1393
1393
!! l_d=(l_{up}l_{down})^{1/2}
1394
1394
!!\f]
1395
1395
!! where \f$l_{up}\f$ and \f$l_{down}\f$ are the distances that a parcel
@@ -1413,7 +1413,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
1413
1413
!
1414
1414
enddo
1415
1415
enddo
1416
- !> - Compute the surface layer length scale (\f$l_1\f$) following
1416
+ !> - Compute the surface layer length scale (\f$l_1\f$) following
1417
1417
!! Nakanishi (2001) \cite Nakanish_2001 (eqn 9 of Han et al.(2019) \cite Han_2019)
1418
1418
do k = 1, km1
1419
1419
do i = 1, im
@@ -1581,7 +1581,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
1581
1581
enddo
1582
1582
enddo
1583
1583
!
1584
- ! eddy diffusivities at model interface (zm level) in LES scale
1584
+ ! eddy diffusivities at model interface (zm level) in LES scale
1585
1585
!
1586
1586
do k = 1, km1
1587
1587
do i = 1, im
@@ -1614,7 +1614,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
1614
1614
do k = 2, kmpbl
1615
1615
do i = 1, im
1616
1616
if(tke(i,k) > tkemax(i)) then
1617
- tkemax(i) = tke(i,k)
1617
+ tkemax(i) = tke(i,k)
1618
1618
ktkemax(i) = k
1619
1619
endif
1620
1620
enddo
@@ -1673,7 +1673,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
1673
1673
enddo
1674
1674
enddo
1675
1675
!
1676
- ! eddy diffusivities at model layer (zl level) in LES scale
1676
+ ! eddy diffusivities at model layer (zl level) in LES scale
1677
1677
!
1678
1678
do k = 1, km1
1679
1679
do i = 1, im
@@ -1717,7 +1717,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
1717
1717
1718
1718
!PCC_CANOPY------------------------------------
1719
1719
kount=0 !IVAI
1720
- if (do_canopy .and. cplaqm) then
1720
+ if (do_canopy .and. cplaqm) then
1721
1721
1722
1722
!IVAI
1723
1723
! print*, ' SATMEDMFVDIFQ_RUN: CLAIE = ' , claie(:)
@@ -1774,7 +1774,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
1774
1774
! & .OR. MAX(0.0, 1.0 - XCANOPYFRT) .GT. 0.5
1775
1775
! & .OR. POPU .GT. 10000.0
1776
1776
! & .OR. EXP(-0.5*XCANOPYLAI*XCANOPYCLU).GT. 0.45
1777
- ! & .AND. FCH .LT. 18.0 ) THEN
1777
+ ! & .AND. FCH .LT. 18.0 ) THEN
1778
1778
1779
1779
dkt(i,k)= dkt(i,k)
1780
1780
dkq(i,k)= dkq(i,k)
@@ -1811,9 +1811,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
1811
1811
IF ( ZCAN/FCH .GT. 1.25 ) THEN !SIGMACAN = Eulerian vertical velocity variance
1812
1812
SIGMACAN = 1.25*ustar(i)
1813
1813
END IF
1814
- IF ( ZCAN/FCH .GE. 0.175
1814
+ IF ( ZCAN/FCH .GE. 0.175
1815
1815
& .AND. ZCAN/FCH .LE. 1.25 ) THEN
1816
- SIGMACAN = ustar(i) * ( 0.75 +
1816
+ SIGMACAN = ustar(i) * ( 0.75 +
1817
1817
& (0.5 * COS((PI/1.06818) *
1818
1818
& (1.25 - (ZCAN/FCH)))) )
1819
1819
END IF
@@ -1825,9 +1825,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
1825
1825
IF ( ZCAN/FCH .GT. 1.25 ) THEN
1826
1826
SIGMACAN = 1.0*ustar(i)
1827
1827
END IF
1828
- IF ( ZCAN/FCH .GE. 0.175
1828
+ IF ( ZCAN/FCH .GE. 0.175
1829
1829
& .AND. ZCAN/FCH .LE. 1.25 ) THEN
1830
- SIGMACAN = ustar(i) * ( 0.625 +
1830
+ SIGMACAN = ustar(i) * ( 0.625 +
1831
1831
& (0.375* COS((PI/1.06818) *
1832
1832
& (1.25 - (ZCAN/FCH)))) )
1833
1833
END IF
@@ -1839,12 +1839,12 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
1839
1839
IF ( ZCAN/FCH .GT. 1.25 ) THEN
1840
1840
SIGMACAN = 0.25*(4.375 - (3.75*HOL))*ustar(i)
1841
1841
END IF
1842
- IF ( ZCAN/FCH .GE. 0.175
1842
+ IF ( ZCAN/FCH .GE. 0.175
1843
1843
& .AND. ZCAN/FCH .LE. 1.25 ) THEN
1844
1844
RRCAN=4.375-(3.75*HOL)
1845
1845
AACAN=(0.125*RRCAN) + 0.125
1846
1846
BBCAN=(0.125*RRCAN) - 0.125
1847
- SIGMACAN = ustar(i) * ( AACAN +
1847
+ SIGMACAN = ustar(i) * ( AACAN +
1848
1848
& (BBCAN * COS((PI/1.06818) *
1849
1849
& (1.25 - (ZCAN/FCH)))) )
1850
1850
END IF
@@ -2051,7 +2051,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
2051
2051
enddo
2052
2052
!
2053
2053
!----------------------------------------------------------------------
2054
- !> - First predict te due to te production & dissipation(diss)
2054
+ !> - First predict te due to te production & dissipation(diss)
2055
2055
!
2056
2056
if(sa3dtke) then
2057
2057
!The following is for SA-3D-TKE
@@ -3052,10 +3052,10 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
3052
3052
enddo
3053
3053
enddo
3054
3054
do i = 1,im
3055
- if(icplocn2atm == 0 ) then
3055
+ if(.not. use_oceanuv ) then
3056
3056
dusfc(i) = -1.*rho_a(i)*stress(i)*u1(i,1)/spd1(i)
3057
3057
dvsfc(i) = -1.*rho_a(i)*stress(i)*v1(i,1)/spd1(i)
3058
- else if (icplocn2atm ==1 ) then
3058
+ else if (use_oceanuv ) then
3059
3059
spd1_m=sqrt( (u1(i,1)-usfco(i))**2+(v1(i,1)-vsfco(i))**2 )
3060
3060
dusfc(i) = -1.*rho_a(i)*stress(i)*(u1(i,1)-usfco(i))/spd1_m
3061
3061
dvsfc(i) = -1.*rho_a(i)*stress(i)*(v1(i,1)-vsfco(i))/spd1_m
0 commit comments