Ignore:
Timestamp:
Nov 21, 2019, 4:43:45 PM (4 years ago)
Author:
lguez
Message:

Merge revisions 3427:3600 of trunk into branch Ocean_skin

Location:
LMDZ6/branches/Ocean_skin
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Ocean_skin

  • LMDZ6/branches/Ocean_skin/libf/phylmd/cloudth_mod.F90

    r2960 r3605  
    77       SUBROUTINE cloudth(ngrid,klev,ind2,  &
    88     &           ztv,po,zqta,fraca, &
    9      &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
     9     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
    1010     &           ratqs,zqs,t)
    1111
    … …  
    3838      REAL zpspsk(ngrid,klev)
    3939      REAL paprs(ngrid,klev+1)
     40      REAL pplay(ngrid,klev)
    4041      REAL ztla(ngrid,klev)
    4142      REAL zthl(ngrid,klev)
    … …  
    7879      CALL cloudth_vert(ngrid,klev,ind2,  &
    7980     &           ztv,po,zqta,fraca, &
    80      &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
     81     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
    8182     &           ratqs,zqs,t)
    8283      RETURN
    … …  
    251252     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
    252253     &           ztv,po,zqta,fraca, &
    253      &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
     254     &           qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
    254255     &           ratqs,zqs,t)
    255256
    … …  
    282283      REAL zpspsk(ngrid,klev)
    283284      REAL paprs(ngrid,klev+1)
     285      REAL pplay(ngrid,klev)
    284286      REAL ztla(ngrid,klev)
    285287      REAL zthl(ngrid,klev)
    … …  
    585587END SUBROUTINE cloudth_vert
    586588
     589
     590
     591
    587592       SUBROUTINE cloudth_v3(ngrid,klev,ind2,  &
    588593     &           ztv,po,zqta,fraca, &
    589      &           qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
     594     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
    590595     &           ratqs,zqs,t)
    591596
    … …  
    618623      REAL zpspsk(ngrid,klev)
    619624      REAL paprs(ngrid,klev+1)
     625      REAL pplay(ngrid,klev)
    620626      REAL ztla(ngrid,klev)
    621627      REAL zthl(ngrid,klev)
    … …  
    641647      REAL alth,alenv,ath,aenv
    642648      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
     649      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
    643650      REAL Tbef,zdelta,qsatbef,zcor
    644651      REAL qlbef 
    … …  
    654661      CALL cloudth_vert_v3(ngrid,klev,ind2,  &
    655662     &           ztv,po,zqta,fraca, &
    656      &           qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
     663     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
    657664     &           ratqs,zqs,t)
    658665      RETURN
    … …  
    808815     SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2,  &
    809816     &           ztv,po,zqta,fraca, &
    810      &           qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, &
     817     &           qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
    811818     &           ratqs,zqs,t)
    812819
    … …  
    841848      REAL zpspsk(ngrid,klev)
    842849      REAL paprs(ngrid,klev+1)
     850      REAL pplay(ngrid,klev)
    843851      REAL ztla(ngrid,klev)
    844852      REAL zthl(ngrid,klev)
    … …  
    864872      REAL alth,alenv,ath,aenv
    865873      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
     874      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
    866875      REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
    867876      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
    … …  
    876885      REAL,SAVE :: sigma1s_factor=1.1
    877886      REAL,SAVE :: sigma1s_power=0.6
     887      REAL,SAVE :: sigma2s_factor=0.09
     888      REAL,SAVE :: sigma2s_power=0.5
    878889      REAL,SAVE :: cloudth_ratqsmin=-1.
    879       !$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,cloudth_ratqsmin)
     890      !$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin)
    880891      INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0
    881892      !$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs)
    … …  
    888899      REAL zqs(ngrid), qcloud(ngrid)
    889900      REAL erf
     901
     902      REAL rhodz(ngrid,klev)
     903      REAL zrho(ngrid,klev)
     904      REAL dz(ngrid,klev)
     905
     906      DO ind1 = 1, ngrid
     907        !Layer calculation
     908        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2
     909        zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3
     910        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre
     911      END DO
     912
    890913
    891914!------------------------------------------------------------------------------
    … …  
    930953        CALL getin_p('cloudth_sigma1s_power',sigma1s_power)
    931954        WRITE(*,*) 'cloudth_sigma1s_power = ', sigma1s_power
     955        ! Factor used in the calculation of sigma2s
     956        CALL getin_p('cloudth_sigma2s_factor',sigma2s_factor)
     957        WRITE(*,*) 'cloudth_sigma2s_factor = ', sigma2s_factor
     958        ! Power used in the calculation of sigma2s
     959        CALL getin_p('cloudth_sigma2s_power',sigma2s_power)
     960        WRITE(*,*) 'cloudth_sigma2s_power = ', sigma2s_power
    932961        ! Minimum value for the environmental air subgrid water distrib
    933962        CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin)
    … …  
    9981027      ENDIF
    9991028      sigma1s = sigma1s_fraca + sigma1s_ratqs
    1000       sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
     1029      sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
    10011030!      tests
    10021031!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
    … …  
    10501079
    10511080      ELSE IF (iflag_cloudth_vert >= 3) THEN
    1052 
     1081      IF (iflag_cloudth_vert < 5) THEN
    10531082!-------------------------------------------------------------------------------
    10541083!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
    … …  
    11191148      endif
    11201149
    1121 
    11221150      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    11231151      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
    11241152
     1153      ELSE IF (iflag_cloudth_vert == 5) THEN
     1154      sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5)+ratqs(ind1,ind2)*po(ind1) !Environment
     1155      sigma2s=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.02)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2)                   !Thermals
     1156      !sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
     1157      !sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
     1158      xth=sth/(sqrt(2.)*sigma2s)
     1159      xenv=senv/(sqrt(2.)*sigma1s)
     1160
     1161      !Volumique
     1162      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
     1163      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
     1164      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
     1165      !print *,'jeanjean_CV=',ctot_vol(ind1,ind2)
     1166
     1167      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth_vol(ind1,ind2))
     1168      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv_vol(ind1,ind2)) 
     1169      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
     1170
     1171      !Surfacique
     1172      !Neggers
     1173      !beta=0.0044
     1174      !inverse_rho=1.+beta*dz(ind1,ind2)
     1175      !print *,'jeanjean : beta=',beta
     1176      !cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
     1177      !cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
     1178      !ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
     1179
     1180      !Brooks
     1181      a_Brooks=0.6694
     1182      b_Brooks=0.1882
     1183      A_Maj_Brooks=0.1635 !-- sans shear
     1184      !A_Maj_Brooks=0.17   !-- ARM LES
     1185      !A_Maj_Brooks=0.18   !-- RICO LES
     1186      !A_Maj_Brooks=0.19   !-- BOMEX LES
     1187      Dx_Brooks=200000.
     1188      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
     1189      !print *,'jeanjean_f=',f_Brooks
     1190
     1191      cth(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
     1192      cenv(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
     1193      ctot(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
     1194      !print *,'JJ_ctot_1',ctot(ind1,ind2)
     1195
     1196
     1197
     1198
     1199
     1200      ENDIF ! of if (iflag_cloudth_vert<5)
    11251201      ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
    11261202
    1127 
     1203!      if (ctot(ind1,ind2).lt.1.e-10) then
    11281204      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
    11291205      ctot(ind1,ind2)=0.
    11301206      ctot_vol(ind1,ind2)=0.
    1131       qcloud(ind1)=zqsatenv(ind1,ind2) 
    1132 
    1133       else 
     1207      qcloud(ind1)=zqsatenv(ind1,ind2)
     1208
     1209      else
    11341210               
    11351211      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
    … …  
    11391215      endif 
    11401216
    1141       else  ! Environment only
     1217      else  ! gaussienne environnement seule
    11421218     
    11431219      zqenv(ind1)=po(ind1)
    … …  
    11511227     
    11521228
    1153 !     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
     1229!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
    11541230      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
    1155       alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
     1231      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
    11561232      aenv=1./(1.+(alenv*Lv/cppd))
    11571233      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
    11581234      sth=0.
     1235     
    11591236
    11601237      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
    11611238      sigma2s=0.
    11621239
    1163       xenv=senv/(sqrt2*sigma1s)
     1240      sqrt2pi=sqrt(2.*pi)
     1241      xenv=senv/(sqrt(2.)*sigma1s)
    11641242      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    11651243      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
    1166       qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
     1244      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
    11671245     
    11681246      if (ctot(ind1,ind2).lt.1.e-3) then
    11691247      ctot(ind1,ind2)=0.
    1170       qcloud(ind1)=zqsatenv(ind1,ind2) 
     1248      qcloud(ind1)=zqsatenv(ind1,ind2)
    11711249
    11721250      else   
    11731251               
     1252!      ctot(ind1,ind2)=ctot(ind1,ind2)
    11741253      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
    11751254
    1176       endif   
    1177  
     1255      endif 
     1256 
     1257
     1258
     1259
    11781260      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
    11791261      ! Outputs used to check the PDFs
    … …  
    11841266
    11851267      enddo       ! from the loop on ngrid l.333
    1186      
    11871268      return
    11881269!     end
    11891270END SUBROUTINE cloudth_vert_v3
    11901271!
     1272
     1273
     1274
     1275
     1276
     1277
     1278
     1279
     1280
     1281
     1282
     1283       SUBROUTINE cloudth_v6(ngrid,klev,ind2,  &
     1284     &           ztv,po,zqta,fraca, &
     1285     &           qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
     1286     &           ratqs,zqs,T)
     1287
     1288
     1289      USE ioipsl_getin_p_mod, ONLY : getin_p
     1290      USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
     1291     &                                cloudth_sigmath,cloudth_sigmaenv
     1292
     1293      IMPLICIT NONE
     1294
     1295#include "YOMCST.h"
     1296#include "YOETHF.h"
     1297#include "FCTTRE.h"
     1298#include "thermcell.h"
     1299#include "nuage.h"
     1300
     1301
     1302        !Domain variables
     1303      INTEGER ngrid !indice Max lat-lon
     1304      INTEGER klev  !indice Max alt
     1305      INTEGER ind1  !indice in [1:ngrid]
     1306      INTEGER ind2  !indice in [1:klev]
     1307        !thermal plume fraction
     1308      REAL fraca(ngrid,klev+1)   !thermal plumes fraction in the gridbox
     1309        !temperatures
     1310      REAL T(ngrid,klev)       !temperature
     1311      REAL zpspsk(ngrid,klev)  !factor (p/p0)**kappa (used for potential variables)
     1312      REAL ztv(ngrid,klev)     !potential temperature (voir thermcell_env.F90)
     1313      REAL ztla(ngrid,klev)    !liquid temperature in the thermals (Tl_th)
     1314      REAL zthl(ngrid,klev)    !liquid temperature in the environment (Tl_env)
     1315        !pressure
     1316      REAL paprs(ngrid,klev+1)   !pressure at the interface of levels
     1317      REAL pplay(ngrid,klev)     !pressure at the middle of the level
     1318        !humidity
     1319      REAL ratqs(ngrid,klev)   !width of the total water subgrid-scale distribution
     1320      REAL po(ngrid)           !total water (qt)
     1321      REAL zqenv(ngrid)        !total water in the environment (qt_env)
     1322      REAL zqta(ngrid,klev)    !total water in the thermals (qt_th)
     1323      REAL zqsatth(ngrid,klev)   !water saturation level in the thermals (q_sat_th)
     1324      REAL zqsatenv(ngrid,klev)  !water saturation level in the environment (q_sat_env)
     1325      REAL qlth(ngrid,klev)    !condensed water in the thermals
     1326      REAL qlenv(ngrid,klev)   !condensed water in the environment
     1327      REAL qltot(ngrid,klev)   !condensed water in the gridbox
     1328        !cloud fractions
     1329      REAL cth_vol(ngrid,klev)   !cloud fraction by volume in the thermals
     1330      REAL cenv_vol(ngrid,klev)  !cloud fraction by volume in the environment
     1331      REAL ctot_vol(ngrid,klev)  !cloud fraction by volume in the gridbox
     1332      REAL cth_surf(ngrid,klev)  !cloud fraction by surface in the thermals
     1333      REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment 
     1334      REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
     1335        !PDF of saturation deficit variables
     1336      REAL rdd,cppd,Lv
     1337      REAL Tbef,zdelta,qsatbef,zcor
     1338      REAL alth,alenv,ath,aenv
     1339      REAL sth,senv              !saturation deficits in the thermals and environment
     1340      REAL sigma_env,sigma_th    !standard deviations of the biGaussian PDF
     1341        !cloud fraction variables
     1342      REAL xth,xenv
     1343      REAL inverse_rho,beta                                  !Neggers et al. (2011) method
     1344      REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
     1345        !Incloud total water variables
     1346      REAL zqs(ngrid)    !q_sat
     1347      REAL qcloud(ngrid) !eau totale dans le nuage
     1348        !Some arithmetic variables
     1349      REAL erf,pi,sqrt2,sqrt2pi
     1350        !Depth of the layer
     1351      REAL dz(ngrid,klev)    !epaisseur de la couche en metre
     1352      REAL rhodz(ngrid,klev)
     1353      REAL zrho(ngrid,klev)
     1354      DO ind1 = 1, ngrid
     1355        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
     1356        zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd          ![kg/m3]
     1357        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2)            ![m]
     1358      END DO
     1359
     1360!------------------------------------------------------------------------------
     1361! Initialization
     1362!------------------------------------------------------------------------------
     1363      qlth(:,:)=0.
     1364      qlenv(:,:)=0. 
     1365      qltot(:,:)=0.
     1366      cth_vol(:,:)=0.
     1367      cenv_vol(:,:)=0.
     1368      ctot_vol(:,:)=0.
     1369      cth_surf(:,:)=0.
     1370      cenv_surf(:,:)=0.
     1371      ctot_surf(:,:)=0.
     1372      qcloud(:)=0.
     1373      rdd=287.04
     1374      cppd=1005.7
     1375      pi=3.14159
     1376      Lv=2.5e6
     1377      sqrt2=sqrt(2.)
     1378      sqrt2pi=sqrt(2.*pi)
     1379
     1380
     1381      DO ind1=1,ngrid
     1382!-------------------------------------------------------------------------------
     1383!Both thermal and environment in the gridbox
     1384!-------------------------------------------------------------------------------
     1385      IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN
     1386        !--------------------------------------------
     1387        !calcul de qsat_env
     1388        !--------------------------------------------
     1389      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
     1390      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
     1391      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     1392      qsatbef=MIN(0.5,qsatbef)
     1393      zcor=1./(1.-retv*qsatbef)
     1394      qsatbef=qsatbef*zcor
     1395      zqsatenv(ind1,ind2)=qsatbef
     1396        !--------------------------------------------
     1397        !calcul de s_env
     1398        !--------------------------------------------
     1399      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
     1400      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
     1401      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam
     1402        !--------------------------------------------
     1403        !calcul de qsat_th
     1404        !--------------------------------------------
     1405      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
     1406      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
     1407      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     1408      qsatbef=MIN(0.5,qsatbef)
     1409      zcor=1./(1.-retv*qsatbef)
     1410      qsatbef=qsatbef*zcor
     1411      zqsatth(ind1,ind2)=qsatbef
     1412        !--------------------------------------------
     1413        !calcul de s_th 
     1414        !--------------------------------------------
     1415      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84 these Arnaud Jam
     1416      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84 these Arnaud Jam
     1417      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84 these Arnaud Jam
     1418        !--------------------------------------------
     1419        !calcul standard deviations bi-Gaussian PDF
     1420        !--------------------------------------------
     1421      sigma_th=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.01)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2)
     1422      sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5)+ratqs(ind1,ind2)*po(ind1)
     1423      xth=sth/(sqrt2*sigma_th)
     1424      xenv=senv/(sqrt2*sigma_env)
     1425        !--------------------------------------------
     1426        !Cloud fraction by volume CF_vol
     1427        !--------------------------------------------
     1428      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
     1429      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
     1430      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
     1431        !--------------------------------------------
     1432        !Condensed water qc
     1433        !--------------------------------------------
     1434      qlth(ind1,ind2)=sigma_th*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth_vol(ind1,ind2))
     1435      qlenv(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv_vol(ind1,ind2)) 
     1436      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
     1437        !--------------------------------------------
     1438        !Cloud fraction by surface CF_surf
     1439        !--------------------------------------------
     1440        !Method Neggers et al. (2011) : ok for cumulus clouds only
     1441      !beta=0.0044 (Jouhaud et al.2018)
     1442      !inverse_rho=1.+beta*dz(ind1,ind2)
     1443      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
     1444        !Method Brooks et al. (2005) : ok for all types of clouds
     1445      a_Brooks=0.6694
     1446      b_Brooks=0.1882
     1447      A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
     1448      Dx_Brooks=200000.   !-- si l'on considere des mailles de 200km de cote
     1449      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
     1450      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
     1451        !--------------------------------------------
     1452        !Incloud Condensed water qcloud
     1453        !--------------------------------------------
     1454      if (ctot_surf(ind1,ind2) .lt. 1.e-10) then
     1455      ctot_vol(ind1,ind2)=0.
     1456      ctot_surf(ind1,ind2)=0.
     1457      qcloud(ind1)=zqsatenv(ind1,ind2)
     1458      else
     1459      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
     1460      endif
     1461
     1462
     1463
     1464!-------------------------------------------------------------------------------
     1465!Environment only in the gridbox
     1466!-------------------------------------------------------------------------------
     1467      ELSE
     1468        !--------------------------------------------
     1469        !calcul de qsat_env
     1470        !--------------------------------------------
     1471      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
     1472      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
     1473      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     1474      qsatbef=MIN(0.5,qsatbef)
     1475      zcor=1./(1.-retv*qsatbef)
     1476      qsatbef=qsatbef*zcor
     1477      zqsatenv(ind1,ind2)=qsatbef
     1478        !--------------------------------------------
     1479        !calcul de s_env
     1480        !--------------------------------------------
     1481      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
     1482      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
     1483      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam
     1484        !--------------------------------------------
     1485        !calcul standard deviations Gaussian PDF
     1486        !--------------------------------------------
     1487      zqenv(ind1)=po(ind1)
     1488      sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
     1489      xenv=senv/(sqrt2*sigma_env)
     1490        !--------------------------------------------
     1491        !Cloud fraction by volume CF_vol
     1492        !--------------------------------------------
     1493      ctot_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
     1494        !--------------------------------------------
     1495        !Condensed water qc
     1496        !--------------------------------------------
     1497      qltot(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*ctot_vol(ind1,ind2))
     1498        !--------------------------------------------
     1499        !Cloud fraction by surface CF_surf
     1500        !--------------------------------------------
     1501        !Method Neggers et al. (2011) : ok for cumulus clouds only
     1502      !beta=0.0044 (Jouhaud et al.2018)
     1503      !inverse_rho=1.+beta*dz(ind1,ind2)
     1504      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
     1505        !Method Brooks et al. (2005) : ok for all types of clouds
     1506      a_Brooks=0.6694
     1507      b_Brooks=0.1882
     1508      A_Maj_Brooks=0.1635 !-- sans dependence au shear
     1509      Dx_Brooks=200000.
     1510      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
     1511      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
     1512        !--------------------------------------------
     1513        !Incloud Condensed water qcloud
     1514        !--------------------------------------------
     1515      if (ctot_surf(ind1,ind2) .lt. 1.e-8) then
     1516      ctot_vol(ind1,ind2)=0.
     1517      ctot_surf(ind1,ind2)=0.
     1518      qcloud(ind1)=zqsatenv(ind1,ind2)
     1519      else
     1520      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
     1521      endif
     1522
     1523
     1524      END IF  ! From the separation (thermal/envrionnement) et (environnement only)
     1525
     1526      ! Outputs used to check the PDFs
     1527      cloudth_senv(ind1,ind2) = senv
     1528      cloudth_sth(ind1,ind2) = sth
     1529      cloudth_sigmaenv(ind1,ind2) = sigma_env
     1530      cloudth_sigmath(ind1,ind2) = sigma_th
     1531
     1532      END DO  ! From the loop on ngrid
     1533      return
     1534
     1535END SUBROUTINE cloudth_v6
    11911536END MODULE cloudth_mod
     1537
     1538
     1539
     1540
Note: See TracChangeset for help on using the changeset viewer.