Ignore:
Timestamp:
Oct 28, 2016, 5:29:26 PM (8 years ago)
Author:
musat
Message:

Oups I missed the reverse update for cloudth; back to svn2662

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing/libf/phylmd/cloudth.F90

    r2669 r2689  
     1!
     2! $Id: cloudth.F90 2689 2016-10-28 17:26:02Z musat $
     3!
    14       SUBROUTINE cloudth(ngrid,klev,ind2,  &
    25     &           ztv,po,zqta,fraca, &
    … …  
    912
    1013!===========================================================================
    11 ! Author : Arnaud Octavio Jam (LMD/CNRS)
     14! Auteur : Arnaud Octavio Jam (LMD/CNRS)
    1215! Date : 25 Mai 2010
    1316! Objet : calcule les valeurs de qc et rneb dans les thermiques
    … …  
    4952      REAL rneb(ngrid,klev)
    5053      REAL t(ngrid,klev)
    51       REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
     54      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
    5255      REAL rdd,cppd,Lv
    5356      REAL alth,alenv,ath,aenv
    54       REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
     57      REAL sth,senv,sigma1s,sigma2s,xth,xenv
    5558      REAL Tbef,zdelta,qsatbef,zcor
    5659      REAL qlbef 
    57       REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
     60      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
     61     
    5862      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
    5963      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
    … …  
    6266
    6367
     68
     69
     70!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     71! Gestion de deux versions de cloudth
     72!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    6473
    6574      IF (iflag_cloudth_vert.GE.1) THEN
    … …  
    7685! Initialisation des variables r?elles
    7786!-------------------------------------------------------------------------------
     87      sigma1(:,:)=0.
     88      sigma2(:,:)=0.
     89      qlth(:,:)=0.
     90      qlenv(:,:)=0. 
     91      qltot(:,:)=0.
     92      rneb(:,:)=0.
     93      qcloud(:)=0.
     94      cth(:,:)=0.
     95      cenv(:,:)=0.
     96      ctot(:,:)=0.
     97      qsatmmussig1=0.
     98      qsatmmussig2=0.
     99      rdd=287.04
     100      cppd=1005.7
     101      pi=3.14159
     102      Lv=2.5e6
     103      sqrt2pi=sqrt(2.*pi)
     104
     105
     106
     107!-------------------------------------------------------------------------------
     108! Calcul de la fraction du thermique et des ?cart-types des distributions
     109!-------------------------------------------------------------------------------                 
     110      do ind1=1,ngrid
     111
     112      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
     113
     114      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
     115
     116
     117!      zqenv(ind1)=po(ind1)
     118      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
     119      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
     120      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     121      qsatbef=MIN(0.5,qsatbef)
     122      zcor=1./(1.-retv*qsatbef)
     123      qsatbef=qsatbef*zcor
     124      zqsatenv(ind1,ind2)=qsatbef
     125
     126
     127
     128
     129      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
     130      aenv=1./(1.+(alenv*Lv/cppd))
     131      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
     132
     133
     134
     135
     136      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
     137      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
     138      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     139      qsatbef=MIN(0.5,qsatbef)
     140      zcor=1./(1.-retv*qsatbef)
     141      qsatbef=qsatbef*zcor
     142      zqsatth(ind1,ind2)=qsatbef
     143           
     144      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
     145      ath=1./(1.+(alth*Lv/cppd))
     146      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
     147     
     148     
     149
     150!------------------------------------------------------------------------------
     151! Calcul des ?cart-types pour s
     152!------------------------------------------------------------------------------
     153
     154!      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
     155!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2)
     156!       if (paprs(ind1,ind2).gt.90000) then
     157!       ratqs(ind1,ind2)=0.002
     158!       else
     159!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
     160!       endif
     161       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
     162       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
     163!       sigma1s=ratqs(ind1,ind2)*po(ind1)
     164!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
     165 
     166!------------------------------------------------------------------------------
     167! Calcul de l'eau condens?e et de la couverture nuageuse
     168!------------------------------------------------------------------------------
     169      sqrt2pi=sqrt(2.*pi)
     170      xth=sth/(sqrt(2.)*sigma2s)
     171      xenv=senv/(sqrt(2.)*sigma1s)
     172      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
     173      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
     174      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
     175
     176      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
     177      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
     178      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
     179
     180!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     181      if (ctot(ind1,ind2).lt.1.e-10) then
     182      ctot(ind1,ind2)=0.
     183      qcloud(ind1)=zqsatenv(ind1,ind2)
     184
     185      else   
     186               
     187      ctot(ind1,ind2)=ctot(ind1,ind2)
     188      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
     189
     190      endif                           
     191     
     192         
     193!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
     194
     195
     196      else  ! gaussienne environnement seule
     197     
     198      zqenv(ind1)=po(ind1)
     199      Tbef=t(ind1,ind2)
     200      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
     201      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     202      qsatbef=MIN(0.5,qsatbef)
     203      zcor=1./(1.-retv*qsatbef)
     204      qsatbef=qsatbef*zcor
     205      zqsatenv(ind1,ind2)=qsatbef
     206     
     207
     208!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
     209      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
     210      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
     211      aenv=1./(1.+(alenv*Lv/cppd))
     212      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
     213     
     214
     215      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
     216
     217      sqrt2pi=sqrt(2.*pi)
     218      xenv=senv/(sqrt(2.)*sigma1s)
     219      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
     220      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
     221     
     222      if (ctot(ind1,ind2).lt.1.e-3) then
     223      ctot(ind1,ind2)=0.
     224      qcloud(ind1)=zqsatenv(ind1,ind2)
     225
     226      else   
     227               
     228      ctot(ind1,ind2)=ctot(ind1,ind2)
     229      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
     230
     231      endif   
     232 
     233 
     234 
     235 
     236 
     237 
     238      endif   
     239      enddo
     240     
     241      return
     242      end
     243
     244
     245
     246!===========================================================================
     247     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
     248     &           ztv,po,zqta,fraca, &
     249     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
     250     &           ratqs,zqs,t)
     251
     252!===========================================================================
     253! Auteur : Arnaud Octavio Jam (LMD/CNRS)
     254! Date : 25 Mai 2010
     255! Objet : calcule les valeurs de qc et rneb dans les thermiques
     256!===========================================================================
     257
     258
     259      USE ioipsl_getin_p_mod, ONLY : getin_p
     260
     261      IMPLICIT NONE
     262
     263#include "YOMCST.h"
     264#include "YOETHF.h"
     265#include "FCTTRE.h"
     266#include "thermcell.h"
     267#include "nuage.h"
     268     
     269      INTEGER itap,ind1,ind2
     270      INTEGER ngrid,klev,klon,l,ig
     271     
     272      REAL ztv(ngrid,klev)
     273      REAL po(ngrid)
     274      REAL zqenv(ngrid)   
     275      REAL zqta(ngrid,klev)
     276         
     277      REAL fraca(ngrid,klev+1)
     278      REAL zpspsk(ngrid,klev)
     279      REAL paprs(ngrid,klev+1)
     280      REAL ztla(ngrid,klev)
     281      REAL zthl(ngrid,klev)
     282
     283      REAL zqsatth(ngrid,klev)
     284      REAL zqsatenv(ngrid,klev)
     285     
     286     
     287      REAL sigma1(ngrid,klev)                                                         
     288      REAL sigma2(ngrid,klev)
     289      REAL qlth(ngrid,klev)
     290      REAL qlenv(ngrid,klev)
     291      REAL qltot(ngrid,klev)
     292      REAL cth(ngrid,klev) 
     293      REAL cenv(ngrid,klev)   
     294      REAL ctot(ngrid,klev)
     295      REAL rneb(ngrid,klev)
     296      REAL t(ngrid,klev)                                                                 
     297      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
     298      REAL rdd,cppd,Lv,sqrt2,sqrtpi
     299      REAL alth,alenv,ath,aenv
     300      REAL sth,senv,sigma1s,sigma2s,xth,xenv
     301      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
     302      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
     303      REAL Tbef,zdelta,qsatbef,zcor
     304      REAL qlbef 
     305      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
     306      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
     307      ! (J Jouhaud, JL Dufresne, JB Madeleine)
     308      REAL,SAVE :: vert_alpha
     309      LOGICAL, SAVE :: firstcall = .TRUE.
     310     
     311      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
     312      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
     313      REAL zqs(ngrid), qcloud(ngrid)
     314      REAL erf
     315
     316!------------------------------------------------------------------------------
     317! Initialisation des variables r?elles
     318!------------------------------------------------------------------------------
    78319      sigma1(:,:)=0.
    79320      sigma2(:,:)=0.
    … …  
    96337      sqrtpi=sqrt(pi)
    97338
    98 
    99 !-------------------------------------------------------------------------------
    100 ! Cloud fraction in the thermals and standard deviation of the PDFs
     339      IF (firstcall) THEN
     340        vert_alpha=0.5
     341        CALL getin_p('cloudth_vert_alpha',vert_alpha)
     342        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
     343        firstcall=.FALSE.
     344      ENDIF
     345
     346!-------------------------------------------------------------------------------
     347! Calcul de la fraction du thermique et des ?cart-types des distributions
    101348!-------------------------------------------------------------------------------                 
    102349      do ind1=1,ngrid
    … …  
    106353      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
    107354
     355
     356!      zqenv(ind1)=po(ind1)
    108357      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
    109358      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    110       qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     359      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
    111360      qsatbef=MIN(0.5,qsatbef)
    112361      zcor=1./(1.-retv*qsatbef)
    … …  
    115364
    116365
    117       alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
    118       aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
    119       senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
    120 
    121 !po = qt de l'environnement ET des thermique
    122 !zqenv = qt environnement
    123 !zqsatenv = qsat environnement
    124 !zthl = Tl environnement
     366
     367
     368      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
     369      aenv=1./(1.+(alenv*Lv/cppd))
     370      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
     371
     372
    125373
    126374
    … …  
    133381      zqsatth(ind1,ind2)=qsatbef
    134382           
    135       alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
    136       ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
    137       sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
    138      
    139 !zqta = qt thermals
    140 !zqsatth = qsat thermals
    141 !ztla = Tl thermals
    142 
    143 !------------------------------------------------------------------------------
    144 ! s standard deviations
    145 !------------------------------------------------------------------------------
    146 
    147 !     tests
    148 !     sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
    149 !     sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
    150 !     sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
    151 !     final option
    152       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
    153       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
    154  
    155 !------------------------------------------------------------------------------
    156 ! Condensed water and cloud cover
    157 !------------------------------------------------------------------------------
    158       xth=sth/(sqrt2*sigma2s)
    159       xenv=senv/(sqrt2*sigma1s)
    160       cth(ind1,ind2)=0.5*(1.+1.*erf(xth))       !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
    161       cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
    162       ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
    163 
    164       qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
    165       qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
    166       qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    167 
    168       if (ctot(ind1,ind2).lt.1.e-10) then
    169       ctot(ind1,ind2)=0.
    170       qcloud(ind1)=zqsatenv(ind1,ind2)
    171       else
    172       qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
    173       endif
    174 
    175       else  ! Environnement only, follow the if l.110
    176      
    177       zqenv(ind1)=po(ind1)
    178       Tbef=t(ind1,ind2)
    179       zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    180       qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
    181       qsatbef=MIN(0.5,qsatbef)
    182       zcor=1./(1.-retv*qsatbef)
    183       qsatbef=qsatbef*zcor
    184       zqsatenv(ind1,ind2)=qsatbef
    185 
    186 !     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
    187       zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
    188       alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
    189       aenv=1./(1.+(alenv*Lv/cppd))
    190       senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
    191      
    192       sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
    193 
    194       xenv=senv/(sqrt2*sigma1s)
    195       ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    196       qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
    197 
    198       if (ctot(ind1,ind2).lt.1.e-3) then
    199       ctot(ind1,ind2)=0.
    200       qcloud(ind1)=zqsatenv(ind1,ind2)
    201       else   
    202       qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
    203       endif
    204 
    205 
    206       endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
    207       enddo       ! from the loop on ngrid l.108
    208       return
    209       end
    210 
    211 
    212 
    213 !===========================================================================
    214      SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
    215      &           ztv,po,zqta,fraca, &
    216      &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
    217      &           ratqs,zqs,t)
    218 
    219 !===========================================================================
    220 ! Auteur : Arnaud Octavio Jam (LMD/CNRS)
    221 ! Date : 25 Mai 2010
    222 ! Objet : calcule les valeurs de qc et rneb dans les thermiques
    223 !===========================================================================
    224 
    225 
    226       USE ioipsl_getin_p_mod, ONLY : getin_p
    227 
    228       IMPLICIT NONE
    229 
    230 #include "YOMCST.h"
    231 #include "YOETHF.h"
    232 #include "FCTTRE.h"
    233 #include "thermcell.h"
    234 #include "nuage.h"
    235      
    236       INTEGER itap,ind1,ind2
    237       INTEGER ngrid,klev,klon,l,ig
    238      
    239       REAL ztv(ngrid,klev)
    240       REAL po(ngrid)
    241       REAL zqenv(ngrid)   
    242       REAL zqta(ngrid,klev)
    243          
    244       REAL fraca(ngrid,klev+1)
    245       REAL zpspsk(ngrid,klev)
    246       REAL paprs(ngrid,klev+1)
    247       REAL ztla(ngrid,klev)
    248       REAL zthl(ngrid,klev)
    249 
    250       REAL zqsatth(ngrid,klev)
    251       REAL zqsatenv(ngrid,klev)
    252      
    253      
    254       REAL sigma1(ngrid,klev)                                                         
    255       REAL sigma2(ngrid,klev)
    256       REAL qlth(ngrid,klev)
    257       REAL qlenv(ngrid,klev)
    258       REAL qltot(ngrid,klev)
    259       REAL cth(ngrid,klev) 
    260       REAL cenv(ngrid,klev)   
    261       REAL ctot(ngrid,klev)
    262       REAL rneb(ngrid,klev)
    263       REAL t(ngrid,klev)                                                                 
    264       REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
    265       REAL rdd,cppd,Lv
    266       REAL alth,alenv,ath,aenv
    267       REAL sth,senv,sigma1s,sigma2s,xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
    268       REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
    269       REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
    270       REAL Tbef,zdelta,qsatbef,zcor
    271       REAL qlbef 
    272       REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
    273       ! Change the width of the PDF used for vertical subgrid scale heterogeneity
    274       ! (J Jouhaud, JL Dufresne, JB Madeleine)
    275       REAL,SAVE :: vert_alpha
    276       LOGICAL, SAVE :: firstcall = .TRUE.
    277 
    278       REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
    279       REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
    280       REAL zqs(ngrid), qcloud(ngrid)
    281       REAL erf
    282 
    283 !------------------------------------------------------------------------------
    284 ! Initialize
    285 !------------------------------------------------------------------------------
    286       sigma1(:,:)=0.
    287       sigma2(:,:)=0.
    288       qlth(:,:)=0.
    289       qlenv(:,:)=0. 
    290       qltot(:,:)=0.
    291       rneb(:,:)=0.
    292       qcloud(:)=0.
    293       cth(:,:)=0.
    294       cenv(:,:)=0.
    295       ctot(:,:)=0.
    296       qsatmmussig1=0.
    297       qsatmmussig2=0.
    298       rdd=287.04
    299       cppd=1005.7
    300       pi=3.14159
    301       Lv=2.5e6
    302       sqrt2pi=sqrt(2.*pi)
    303       sqrt2=sqrt(2.)
    304       sqrtpi=sqrt(pi)
    305 
    306       IF (firstcall) THEN
    307         vert_alpha=0.5
    308         CALL getin_p('cloudth_vert_alpha',vert_alpha)
    309         WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
    310         firstcall=.FALSE.
    311       ENDIF
    312 
    313 !-------------------------------------------------------------------------------
    314 ! Calcul de la fraction du thermique et des ecart-types des distributions
    315 !-------------------------------------------------------------------------------                 
    316       do ind1=1,ngrid
    317 
    318       if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
    319 
    320       zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
    321 
    322 
    323       Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
    324       zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    325       qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
    326       qsatbef=MIN(0.5,qsatbef)
    327       zcor=1./(1.-retv*qsatbef)
    328       qsatbef=qsatbef*zcor
    329       zqsatenv(ind1,ind2)=qsatbef
    330 
    331 
    332       alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
    333       aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
    334       senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
    335 
    336 !zqenv = qt environnement
    337 !zqsatenv = qsat environnement
    338 !zthl = Tl environnement
    339 
    340 
    341       Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
    342       zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    343       qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
    344       qsatbef=MIN(0.5,qsatbef)
    345       zcor=1./(1.-retv*qsatbef)
    346       qsatbef=qsatbef*zcor
    347       zqsatth(ind1,ind2)=qsatbef
    348            
    349       alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
    350       ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
    351       sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
    352      
    353      
    354 !zqta = qt thermals
    355 !zqsatth = qsat thermals
    356 !ztla = Tl thermals
    357 
    358 !------------------------------------------------------------------------------
    359 ! s standard deviation
    360 !------------------------------------------------------------------------------
    361 
    362       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
    363       sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
    364 !      tests
    365 !      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
    366 !      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
    367 !      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
    368 !      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
     383      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
     384      ath=1./(1.+(alth*Lv/cppd))
     385      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
     386     
     387     
     388
     389!------------------------------------------------------------------------------
     390! Calcul des ?cart-types pour s
     391!------------------------------------------------------------------------------
     392
     393      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
     394      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
    369395!       if (paprs(ind1,ind2).gt.90000) then
    370396!       ratqs(ind1,ind2)=0.002
    … …  
    377403!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
    378404 
     405!------------------------------------------------------------------------------
     406! Calcul de l'eau condens?e et de la couverture nuageuse
     407!------------------------------------------------------------------------------
     408      sqrt2pi=sqrt(2.*pi)
     409      xth=sth/(sqrt(2.)*sigma2s)
     410      xenv=senv/(sqrt(2.)*sigma1s)
     411      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
     412      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
     413      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
     414
     415      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
     416      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
     417      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
     418     
    379419       IF (iflag_cloudth_vert == 1) THEN
    380420!-------------------------------------------------------------------------------
    381 !  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
    382 !-------------------------------------------------------------------------------
    383 
     421!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
     422!-------------------------------------------------------------------------------
     423!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
     424!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
    384425      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
    385426      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
    386 
     427!      deltasenv=aenv*0.01*po(ind1)
     428!     deltasth=ath*0.01*zqta(ind1,ind2)   
    387429      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
    388430      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
    … …  
    396438      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
    397439
    398       ! Environment
    399440      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
    400441      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
    … …  
    403444
    404445      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
    405 
    406       ! Thermal
     446!      qlenv(ind1,ind2)=IntJ
     447!      print*, qlenv(ind1,ind2),'VERIF EAU'
     448
     449
    407450      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
     451!      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
     452!      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
    408453      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
    409454      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
    410455      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
    411456      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
     457!      qlth(ind1,ind2)=IntJ
     458!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
    412459      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    413460
    … …  
    415462
    416463!-------------------------------------------------------------------------------
    417 !  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
     464!  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
    418465!-------------------------------------------------------------------------------
    419466!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
    … …  
    426473      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
    427474      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
    428       exp_xenv1 = exp(-1.*xenv1**2)
    429       exp_xenv2 = exp(-1.*xenv2**2)
    430475      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
    431476      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
    432       exp_xth1 = exp(-1.*xth1**2)
    433       exp_xth2 = exp(-1.*xth2**2)
     477!     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
     478!     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
    434479     
    435480      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
    … …  
    437482      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
    438483
    439      
    440       !environnement
    441       IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
    442       if (deltasenv .lt. 1.e-10) then
    443       qlenv(ind1,ind2)=IntJ
    444       else
     484      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
    445485      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
    446       IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
    447       IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
     486      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
     487      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))
     488
     489!      IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
     490!      IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
     491!      IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
     492
    448493      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
    449       endif
    450 
    451 
    452       !thermique
    453       IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
    454       if (deltasth .lt. 1.e-10) then ! Seuil a choisir !!!
    455       qlth(ind1,ind2)=IntJ
    456       else
     494!      qlenv(ind1,ind2)=IntJ
     495!      print*, qlenv(ind1,ind2),'VERIF EAU'
     496
     497      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
    457498      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
    458       IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
    459       IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
     499      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
     500      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
     501     
    460502      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
    461       endif
    462 
     503!      qlth(ind1,ind2)=IntJ
     504!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
    463505      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
     506     
     507
    464508
    465509
    466510      ENDIF ! of if (iflag_cloudth_vert==1 or 2)
     511
     512!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    467513
    468514      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
    … …  
    472518      else
    473519               
     520      ctot(ind1,ind2)=ctot(ind1,ind2)
    474521      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
    475522!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
    … …  
    477524
    478525      endif 
    479 
    480       else  ! Environment only
     526                       
     527     
     528         
     529!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
     530
     531
     532      else  ! gaussienne environnement seule
    481533     
    482534      zqenv(ind1)=po(ind1)
    … …  
    490542     
    491543
    492 !     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
     544!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
    493545      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
    494546      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
    … …  
    499551      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
    500552
    501       xenv=senv/(sqrt2*sigma1s)
     553      sqrt2pi=sqrt(2.*pi)
     554      xenv=senv/(sqrt(2.)*sigma1s)
    502555      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    503       qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
     556      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
    504557     
    505558      if (ctot(ind1,ind2).lt.1.e-3) then
    … …  
    509562      else   
    510563               
     564      ctot(ind1,ind2)=ctot(ind1,ind2)
    511565      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
    512566
    513567      endif   
    514568 
    515       endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
    516       enddo       ! from the loop on ngrid l.333
     569 
     570 
     571 
     572 
     573 
     574      endif   
     575      enddo
    517576     
    518577      return
Note: See TracChangeset for help on using the changeset viewer.