Ignore:
Timestamp:
Jul 18, 2016, 9:41:10 PM (8 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes r2545:2589 into testing branch

Location:
LMDZ5/branches/testing
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • LMDZ5/branches/testing

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

    r2544 r2594  
    55
    66
    7       USE IOIPSL, ONLY : getin
    87      IMPLICIT NONE
    98
     
    2019#include "FCTTRE.h"
    2120#include "thermcell.h"
     21#include "nuage.h"
    2222
    2323      INTEGER itap,ind1,ind2
     
    5454      REAL sth,senv,sigma1s,sigma2s,xth,xenv
    5555      REAL Tbef,zdelta,qsatbef,zcor
    56       REAL alpha,qlbef 
     56      REAL qlbef 
    5757      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
    5858     
     
    6262      REAL erf
    6363
    64       REAL, SAVE :: iflag_cloudth_vert, iflag_cloudth_vert_omp=0
    65 
    66 
    67       LOGICAL, SAVE :: first=.true.
    68 
    69 
    7064
    7165
    7266
    7367!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    74 ! Astuce pour gérer deux versions de cloudth en attendant
    75 ! de converger sur une version nouvelle
     68! Gestion de deux versions de cloudth
    7669!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    77       IF (first) THEN
    78      !$OMP MASTER
    79      CALL getin('iflag_cloudth_vert',iflag_cloudth_vert_omp)
    80      !$OMP END MASTER
    81      !$OMP BARRIER
    82      iflag_cloudth_vert=iflag_cloudth_vert_omp
    83       first=.false.
    84       ENDIF
    85        IF (iflag_cloudth_vert==1) THEN
    86        CALL cloudth_vert(ngrid,klev,ind2,  &
     70
     71      IF (iflag_cloudth_vert.GE.1) THEN
     72      CALL cloudth_vert(ngrid,klev,ind2,  &
    8773     &           ztv,po,zqta,fraca, &
    8874     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
    8975     &           ratqs,zqs,t)
    90        RETURN
    91        ENDIF
     76      RETURN
     77      ENDIF
    9278!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    93 
    9479
    9580
     
    185170      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    186171      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
    187 !      ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*cenv(ind1,ind2)
    188 
    189 
    190172
    191173      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
    192174      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
    193175      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    194 !      qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2)
    195      
    196 
    197 !      print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
    198176
    199177!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
    269247     &           ratqs,zqs,t)
    270248
    271 
    272       IMPLICIT NONE
    273 
    274 
    275249!===========================================================================
    276250! Auteur : Arnaud Octavio Jam (LMD/CNRS)
     
    280254
    281255
     256      USE ioipsl_getin_p_mod, ONLY : getin_p
     257
     258      IMPLICIT NONE
     259
    282260#include "YOMCST.h"
    283261#include "YOETHF.h"
    284262#include "FCTTRE.h"
    285263#include "thermcell.h"
    286 
     264#include "nuage.h"
     265     
    287266      INTEGER itap,ind1,ind2
    288267      INTEGER ngrid,klev,klon,l,ig
     
    320299      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
    321300      REAL Tbef,zdelta,qsatbef,zcor
    322       REAL alpha,qlbef 
     301      REAL qlbef 
    323302      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.
    324307     
    325308      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
     
    327310      REAL zqs(ngrid), qcloud(ngrid)
    328311      REAL erf
    329 
    330 
    331 
    332 
    333312
    334313!------------------------------------------------------------------------------
     
    355334      sqrtpi=sqrt(pi)
    356335
    357 
     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
    358342
    359343!-------------------------------------------------------------------------------
     
    425409      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    426410      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
    427 !      ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*cenv(ind1,ind2)
    428 
    429 
    430411
    431412      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
    432413      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
    433414      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    434 !      qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2)
    435415     
    436 
    437 !      print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
    438 
    439 
     416       IF (iflag_cloudth_vert == 1) THEN
    440417!-------------------------------------------------------------------------------
    441418!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
     
    479456      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    480457
     458      ELSE IF (iflag_cloudth_vert == 2) THEN
     459
     460!-------------------------------------------------------------------------------
     461!  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
     462!-------------------------------------------------------------------------------
     463!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
     464!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
     465!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
     466!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
     467      deltasenv=aenv*vert_alpha*sigma1s
     468      deltasth=ath*vert_alpha*sigma2s
     469     
     470      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
     471      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
     472      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
     473      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
     474!     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
     475!     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
     476     
     477      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
     478      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
     479      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
     480
     481      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
     482      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
     490      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)
     495      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     
     499      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
     500!      qlth(ind1,ind2)=IntJ
     501!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
     502      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
     503     
     504
     505
     506
     507      ENDIF ! of if (iflag_cloudth_vert==1 or 2)
     508
    481509!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     510
    482511      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
    483512      ctot(ind1,ind2)=0.
Note: See TracChangeset for help on using the changeset viewer.