Changeset 2662


Ignore:
Timestamp:
Oct 11, 2016, 7:23:57 PM (8 years ago)
Author:
jbmadeleine
Message:

Updated cloudth.F90 following changes by J. Jouhaud

  • improved calculation of the PDF for vertical subgrid scale heterogeneity

(used when iflag_cloudth > 0)

  • use of ratqs in the computation of sigma1s (even when iflag_cloudth.EQ.0)

Sigma1s should have depended on ratqs before; this was probably a remnant of previous tests.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/cloudth.F90

    r2586 r2662  
    99
    1010!===========================================================================
    11 ! Auteur : Arnaud Octavio Jam (LMD/CNRS)
     11! Author : Arnaud Octavio Jam (LMD/CNRS)
    1212! Date : 25 Mai 2010
    1313! Objet : calcule les valeurs de qc et rneb dans les thermiques
    … …  
    4949      REAL rneb(ngrid,klev)
    5050      REAL t(ngrid,klev)
    51       REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
     51      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
    5252      REAL rdd,cppd,Lv
    5353      REAL alth,alenv,ath,aenv
    54       REAL sth,senv,sigma1s,sigma2s,xth,xenv
     54      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
    5555      REAL Tbef,zdelta,qsatbef,zcor
    5656      REAL qlbef 
    57       REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
    58      
     57      REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
    5958      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
    6059      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
    … …  
    6362
    6463
    65 
    66 
    67 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    68 ! Gestion de deux versions de cloudth
    69 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    7064
    7165      IF (iflag_cloudth_vert.GE.1) THEN
    … …  
    8276! Initialisation des variables r?elles
    8377!-------------------------------------------------------------------------------
    84       sigma1(:,:)=0.
    85       sigma2(:,:)=0.
    86       qlth(:,:)=0.
    87       qlenv(:,:)=0. 
    88       qltot(:,:)=0.
    89       rneb(:,:)=0.
    90       qcloud(:)=0.
    91       cth(:,:)=0.
    92       cenv(:,:)=0.
    93       ctot(:,:)=0.
    94       qsatmmussig1=0.
    95       qsatmmussig2=0.
    96       rdd=287.04
    97       cppd=1005.7
    98       pi=3.14159
    99       Lv=2.5e6
    100       sqrt2pi=sqrt(2.*pi)
    101 
    102 
    103 
    104 !-------------------------------------------------------------------------------
    105 ! Calcul de la fraction du thermique et des ?cart-types des distributions
    106 !-------------------------------------------------------------------------------                 
    107       do ind1=1,ngrid
    108 
    109       if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
    110 
    111       zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
    112 
    113 
    114 !      zqenv(ind1)=po(ind1)
    115       Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
    116       zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    117       qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
    118       qsatbef=MIN(0.5,qsatbef)
    119       zcor=1./(1.-retv*qsatbef)
    120       qsatbef=qsatbef*zcor
    121       zqsatenv(ind1,ind2)=qsatbef
    122 
    123 
    124 
    125 
    126       alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
    127       aenv=1./(1.+(alenv*Lv/cppd))
    128       senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
    129 
    130 
    131 
    132 
    133       Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
    134       zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    135       qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
    136       qsatbef=MIN(0.5,qsatbef)
    137       zcor=1./(1.-retv*qsatbef)
    138       qsatbef=qsatbef*zcor
    139       zqsatth(ind1,ind2)=qsatbef
    140            
    141       alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
    142       ath=1./(1.+(alth*Lv/cppd))
    143       sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
    144      
    145      
    146 
    147 !------------------------------------------------------------------------------
    148 ! Calcul des ?cart-types pour s
    149 !------------------------------------------------------------------------------
    150 
    151 !      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
    152 !      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2)
    153 !       if (paprs(ind1,ind2).gt.90000) then
    154 !       ratqs(ind1,ind2)=0.002
    155 !       else
    156 !       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
    157 !       endif
    158        sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
    159        sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
    160 !       sigma1s=ratqs(ind1,ind2)*po(ind1)
    161 !      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
    162  
    163 !------------------------------------------------------------------------------
    164 ! Calcul de l'eau condens?e et de la couverture nuageuse
    165 !------------------------------------------------------------------------------
    166       sqrt2pi=sqrt(2.*pi)
    167       xth=sth/(sqrt(2.)*sigma2s)
    168       xenv=senv/(sqrt(2.)*sigma1s)
    169       cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
    170       cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    171       ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
    172 
    173       qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
    174       qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
    175       qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    176 
    177 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    178       if (ctot(ind1,ind2).lt.1.e-10) then
    179       ctot(ind1,ind2)=0.
    180       qcloud(ind1)=zqsatenv(ind1,ind2)
    181 
    182       else   
    183                
    184       ctot(ind1,ind2)=ctot(ind1,ind2)
    185       qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
    186 
    187       endif                           
    188      
    189          
    190 !     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
    191 
    192 
    193       else  ! gaussienne environnement seule
    194      
    195       zqenv(ind1)=po(ind1)
    196       Tbef=t(ind1,ind2)
    197       zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    198       qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
    199       qsatbef=MIN(0.5,qsatbef)
    200       zcor=1./(1.-retv*qsatbef)
    201       qsatbef=qsatbef*zcor
    202       zqsatenv(ind1,ind2)=qsatbef
    203      
    204 
    205 !      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
    206       zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
    207       alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
    208       aenv=1./(1.+(alenv*Lv/cppd))
    209       senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
    210      
    211 
    212       sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
    213 
    214       sqrt2pi=sqrt(2.*pi)
    215       xenv=senv/(sqrt(2.)*sigma1s)
    216       ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    217       qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
    218      
    219       if (ctot(ind1,ind2).lt.1.e-3) then
    220       ctot(ind1,ind2)=0.
    221       qcloud(ind1)=zqsatenv(ind1,ind2)
    222 
    223       else   
    224                
    225       ctot(ind1,ind2)=ctot(ind1,ind2)
    226       qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
    227 
    228       endif   
    229  
    230  
    231  
    232  
    233  
    234  
    235       endif   
    236       enddo
    237      
    238       return
    239       end
    240 
    241 
    242 
    243 !===========================================================================
    244      SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
    245      &           ztv,po,zqta,fraca, &
    246      &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
    247      &           ratqs,zqs,t)
    248 
    249 !===========================================================================
    250 ! Auteur : Arnaud Octavio Jam (LMD/CNRS)
    251 ! Date : 25 Mai 2010
    252 ! Objet : calcule les valeurs de qc et rneb dans les thermiques
    253 !===========================================================================
    254 
    255 
    256       USE ioipsl_getin_p_mod, ONLY : getin_p
    257 
    258       IMPLICIT NONE
    259 
    260 #include "YOMCST.h"
    261 #include "YOETHF.h"
    262 #include "FCTTRE.h"
    263 #include "thermcell.h"
    264 #include "nuage.h"
    265      
    266       INTEGER itap,ind1,ind2
    267       INTEGER ngrid,klev,klon,l,ig
    268      
    269       REAL ztv(ngrid,klev)
    270       REAL po(ngrid)
    271       REAL zqenv(ngrid)   
    272       REAL zqta(ngrid,klev)
    273          
    274       REAL fraca(ngrid,klev+1)
    275       REAL zpspsk(ngrid,klev)
    276       REAL paprs(ngrid,klev+1)
    277       REAL ztla(ngrid,klev)
    278       REAL zthl(ngrid,klev)
    279 
    280       REAL zqsatth(ngrid,klev)
    281       REAL zqsatenv(ngrid,klev)
    282      
    283      
    284       REAL sigma1(ngrid,klev)                                                         
    285       REAL sigma2(ngrid,klev)
    286       REAL qlth(ngrid,klev)
    287       REAL qlenv(ngrid,klev)
    288       REAL qltot(ngrid,klev)
    289       REAL cth(ngrid,klev) 
    290       REAL cenv(ngrid,klev)   
    291       REAL ctot(ngrid,klev)
    292       REAL rneb(ngrid,klev)
    293       REAL t(ngrid,klev)                                                                 
    294       REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
    295       REAL rdd,cppd,Lv,sqrt2,sqrtpi
    296       REAL alth,alenv,ath,aenv
    297       REAL sth,senv,sigma1s,sigma2s,xth,xenv
    298       REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
    299       REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
    300       REAL Tbef,zdelta,qsatbef,zcor
    301       REAL qlbef 
    302       REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
    303       ! Change the width of the PDF used for vertical subgrid scale heterogeneity
    304       ! (J Jouhaud, JL Dufresne, JB Madeleine)
    305       REAL,SAVE :: vert_alpha
    306       LOGICAL, SAVE :: firstcall = .TRUE.
    307      
    308       REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
    309       REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
    310       REAL zqs(ngrid), qcloud(ngrid)
    311       REAL erf
    312 
    313 !------------------------------------------------------------------------------
    314 ! Initialisation des variables r?elles
    315 !------------------------------------------------------------------------------
    31678      sigma1(:,:)=0.
    31779      sigma2(:,:)=0.
    … …  
    33496      sqrtpi=sqrt(pi)
    33597
    336       IF (firstcall) THEN
    337         vert_alpha=0.5
    338         CALL getin_p('cloudth_vert_alpha',vert_alpha)
    339         WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
    340         firstcall=.FALSE.
    341       ENDIF
    342 
    343 !-------------------------------------------------------------------------------
    344 ! Calcul de la fraction du thermique et des ?cart-types des distributions
     98
     99!-------------------------------------------------------------------------------
     100! Cloud fraction in the thermals and standard deviation of the PDFs
    345101!-------------------------------------------------------------------------------                 
    346102      do ind1=1,ngrid
    … …  
    350106      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
    351107
    352 
    353 !      zqenv(ind1)=po(ind1)
    354108      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
    355109      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    356       qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     110      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
    357111      qsatbef=MIN(0.5,qsatbef)
    358112      zcor=1./(1.-retv*qsatbef)
    … …  
    361115
    362116
    363 
    364 
    365       alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
    366       aenv=1./(1.+(alenv*Lv/cppd))
    367       senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
    368 
    369 
     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
    370125
    371126
    … …  
    378133      zqsatth(ind1,ind2)=qsatbef
    379134           
    380       alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
    381       ath=1./(1.+(alth*Lv/cppd))
    382       sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
    383      
    384      
    385 
    386 !------------------------------------------------------------------------------
    387 ! Calcul des ?cart-types pour s
    388 !------------------------------------------------------------------------------
    389 
    390       sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
    391       sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
     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)
    392369!       if (paprs(ind1,ind2).gt.90000) then
    393370!       ratqs(ind1,ind2)=0.002
    … …  
    400377!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
    401378 
    402 !------------------------------------------------------------------------------
    403 ! Calcul de l'eau condens?e et de la couverture nuageuse
    404 !------------------------------------------------------------------------------
    405       sqrt2pi=sqrt(2.*pi)
    406       xth=sth/(sqrt(2.)*sigma2s)
    407       xenv=senv/(sqrt(2.)*sigma1s)
    408       cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
    409       cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    410       ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
    411 
    412       qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
    413       qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
    414       qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    415      
    416379       IF (iflag_cloudth_vert == 1) THEN
    417380!-------------------------------------------------------------------------------
    418 !  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
    419 !-------------------------------------------------------------------------------
    420 !      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
    421 !      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
     381!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
     382!-------------------------------------------------------------------------------
     383
    422384      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
    423385      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
    424 !      deltasenv=aenv*0.01*po(ind1)
    425 !     deltasth=ath*0.01*zqta(ind1,ind2)   
     386
    426387      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
    427388      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
    … …  
    435396      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
    436397
     398      ! Environment
    437399      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
    438400      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
    … …  
    441403
    442404      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
    443 !      qlenv(ind1,ind2)=IntJ
    444 !      print*, qlenv(ind1,ind2),'VERIF EAU'
    445 
    446 
     405
     406      ! Thermal
    447407      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
    448 !      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
    449 !      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
    450408      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
    451409      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
    452410      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
    453411      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
    454 !      qlth(ind1,ind2)=IntJ
    455 !      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
    456412      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    457413
    … …  
    459415
    460416!-------------------------------------------------------------------------------
    461 !  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
     417!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
    462418!-------------------------------------------------------------------------------
    463419!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
    … …  
    470426      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
    471427      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
     428      exp_xenv1 = exp(-1.*xenv1**2)
     429      exp_xenv2 = exp(-1.*xenv2**2)
    472430      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
    473431      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
    474 !     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
    475 !     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
     432      exp_xth1 = exp(-1.*xth1**2)
     433      exp_xth2 = exp(-1.*xth2**2)
    476434     
    477435      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
    … …  
    479437      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
    480438
    481       IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
     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
    482445      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
    483       IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
    484       IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))
    485 
    486 !      IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
    487 !      IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
    488 !      IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
    489 
     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)
    490448      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
    491 !      qlenv(ind1,ind2)=IntJ
    492 !      print*, qlenv(ind1,ind2),'VERIF EAU'
    493 
    494       IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
     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
    495457      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
    496       IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
    497       IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
    498      
     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)
    499460      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
    500 !      qlth(ind1,ind2)=IntJ
    501 !      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
     461      endif
     462
    502463      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    503      
    504 
    505464
    506465
    507466      ENDIF ! of if (iflag_cloudth_vert==1 or 2)
    508 
    509 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    510467
    511468      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
    … …  
    515472      else
    516473               
    517       ctot(ind1,ind2)=ctot(ind1,ind2)
    518474      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
    519475!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
    … …  
    521477
    522478      endif 
    523                        
    524      
    525          
    526 !     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
    527 
    528 
    529       else  ! gaussienne environnement seule
     479
     480      else  ! Environment only
    530481     
    531482      zqenv(ind1)=po(ind1)
    … …  
    539490     
    540491
    541 !      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
     492!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
    542493      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
    543494      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
    … …  
    548499      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
    549500
    550       sqrt2pi=sqrt(2.*pi)
    551       xenv=senv/(sqrt(2.)*sigma1s)
     501      xenv=senv/(sqrt2*sigma1s)
    552502      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    553       qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
     503      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
    554504     
    555505      if (ctot(ind1,ind2).lt.1.e-3) then
    … …  
    559509      else   
    560510               
    561       ctot(ind1,ind2)=ctot(ind1,ind2)
    562511      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
    563512
    564513      endif   
    565514 
    566  
    567  
    568  
    569  
    570  
    571       endif   
    572       enddo
     515      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
     516      enddo       ! from the loop on ngrid l.333
    573517     
    574518      return
Note: See TracChangeset for help on using the changeset viewer.