Ignore:
Timestamp:
Jun 14, 2015, 9:13:32 PM (9 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes -r2237:2291 into testing branch

Location:
LMDZ5/branches/testing
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

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

    r2160 r2298  
    55
    66
     7      USE IOIPSL, ONLY : getin
    78      IMPLICIT NONE
    89
    … …  
    3940     
    4041     
    41       REAL sigma1(ngrid,klev)                                                         
     42      REAL sigma1(ngrid,klev)
    4243      REAL sigma2(ngrid,klev)
    4344      REAL qlth(ngrid,klev)
    … …  
    4849      REAL ctot(ngrid,klev)
    4950      REAL rneb(ngrid,klev)
    50       REAL t(ngrid,klev)                                                                 
     51      REAL t(ngrid,klev)
    5152      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
    5253      REAL rdd,cppd,Lv
    … …  
    6263      REAL erf
    6364
    64 
    65 
    66 
    67 
    68 !      print*,ngrid,klev,ind1,ind2,ztv(ind1,ind2),po(ind1),zqta(ind1,ind2), &
    69 !     &       fraca(ind1,ind2),zpspsk(ind1,ind2),paprs(ind1,ind2),ztla(ind1,ind2),zthl(ind1,ind2), &
    70 !     &       'verif'
    71 
    72 
    73 !      LOGICAL active(ngrid)   
    74      
    75 !-----------------------------------------------------------------------------------------------------------------
     65      REAL, SAVE :: iflag_cloudth_vert, iflag_cloudth_vert_omp=0
     66
     67
     68      LOGICAL, SAVE :: first=.true.
     69
     70
     71
     72
     73
     74!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     75! Astuce pour gérer deux versions de cloudth en attendant
     76! de converger sur une version nouvelle
     77!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     78      IF (first) THEN
     79     !$OMP MASTER
     80     CALL getin('iflag_cloudth_vert',iflag_cloudth_vert_omp)
     81     !$OMP END MASTER
     82     !$OMP BARRIER
     83     iflag_cloudth_vert=iflag_cloudth_vert_omp
     84      first=.false.
     85      ENDIF
     86       IF (iflag_cloudth_vert==1) THEN
     87       CALL cloudth_vert(ngrid,klev,ind2,  &
     88     &           ztv,po,zqta,fraca, &
     89     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
     90     &           ratqs,zqs,t)
     91       RETURN
     92       ENDIF
     93!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     94
     95
     96
     97!-------------------------------------------------------------------------------
    7698! Initialisation des variables r?elles
    77 !-----------------------------------------------------------------------------------------------------------------
     99!-------------------------------------------------------------------------------
    78100      sigma1(:,:)=0.
    79101      sigma2(:,:)=0.
    … …  
    96118
    97119
    98 !------------------------------------------------------------------------------------------------------------------
     120!-------------------------------------------------------------------------------
    99121! Calcul de la fraction du thermique et des ?cart-types des distributions
    100 !------------------------------------------------------------------------------------------------------------------                 
     122!-------------------------------------------------------------------------------                 
    101123      do ind1=1,ngrid
    102124
    … …  
    139161     
    140162
    141 !-----------------------------------------------------------------------------------------------------------------
     163!------------------------------------------------------------------------------
    142164! Calcul des ?cart-types pour s
    143 !-----------------------------------------------------------------------------------------------------------------
     165!------------------------------------------------------------------------------
    144166
    145167!      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
    … …  
    155177!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
    156178 
    157 !-----------------------------------------------------------------------------------------------------------------
     179!------------------------------------------------------------------------------
    158180! Calcul de l'eau condens?e et de la couverture nuageuse
    159 !-----------------------------------------------------------------------------------------------------------------
     181!------------------------------------------------------------------------------
    160182      sqrt2pi=sqrt(2.*pi)
    161183      xth=sth/(sqrt(2.)*sigma2s)
    … …  
    176198!      print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
    177199
    178 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     200!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    179201      if (ctot(ind1,ind2).lt.1.e-10) then
    180202      ctot(ind1,ind2)=0.
    … …  
    242264
    243265
    244 
    245 
    246                                                                            
    247 
    248 
    249 
    250 
     266!===========================================================================
     267     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
     268     &           ztv,po,zqta,fraca, &
     269     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
     270     &           ratqs,zqs,t)
     271
     272
     273      IMPLICIT NONE
     274
     275
     276!===========================================================================
     277! Auteur : Arnaud Octavio Jam (LMD/CNRS)
     278! Date : 25 Mai 2010
     279! Objet : calcule les valeurs de qc et rneb dans les thermiques
     280!===========================================================================
     281
     282
     283#include "YOMCST.h"
     284#include "YOETHF.h"
     285#include "FCTTRE.h"
     286#include "iniprint.h"
     287#include "thermcell.h"
     288
     289      INTEGER itap,ind1,ind2
     290      INTEGER ngrid,klev,klon,l,ig
     291     
     292      REAL ztv(ngrid,klev)
     293      REAL po(ngrid)
     294      REAL zqenv(ngrid)   
     295      REAL zqta(ngrid,klev)
     296         
     297      REAL fraca(ngrid,klev+1)
     298      REAL zpspsk(ngrid,klev)
     299      REAL paprs(ngrid,klev+1)
     300      REAL ztla(ngrid,klev)
     301      REAL zthl(ngrid,klev)
     302
     303      REAL zqsatth(ngrid,klev)
     304      REAL zqsatenv(ngrid,klev)
     305     
     306     
     307      REAL sigma1(ngrid,klev)                                                         
     308      REAL sigma2(ngrid,klev)
     309      REAL qlth(ngrid,klev)
     310      REAL qlenv(ngrid,klev)
     311      REAL qltot(ngrid,klev)
     312      REAL cth(ngrid,klev) 
     313      REAL cenv(ngrid,klev)   
     314      REAL ctot(ngrid,klev)
     315      REAL rneb(ngrid,klev)
     316      REAL t(ngrid,klev)                                                                 
     317      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
     318      REAL rdd,cppd,Lv,sqrt2,sqrtpi
     319      REAL alth,alenv,ath,aenv
     320      REAL sth,senv,sigma1s,sigma2s,xth,xenv
     321      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
     322      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
     323      REAL Tbef,zdelta,qsatbef,zcor
     324      REAL alpha,qlbef 
     325      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
     326     
     327      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
     328      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
     329      REAL zqs(ngrid), qcloud(ngrid)
     330      REAL erf
     331
     332
     333
     334
     335
     336!------------------------------------------------------------------------------
     337! Initialisation des variables r?elles
     338!------------------------------------------------------------------------------
     339      sigma1(:,:)=0.
     340      sigma2(:,:)=0.
     341      qlth(:,:)=0.
     342      qlenv(:,:)=0. 
     343      qltot(:,:)=0.
     344      rneb(:,:)=0.
     345      qcloud(:)=0.
     346      cth(:,:)=0.
     347      cenv(:,:)=0.
     348      ctot(:,:)=0.
     349      qsatmmussig1=0.
     350      qsatmmussig2=0.
     351      rdd=287.04
     352      cppd=1005.7
     353      pi=3.14159
     354      Lv=2.5e6
     355      sqrt2pi=sqrt(2.*pi)
     356      sqrt2=sqrt(2.)
     357      sqrtpi=sqrt(pi)
     358
     359
     360
     361!-------------------------------------------------------------------------------
     362! Calcul de la fraction du thermique et des ?cart-types des distributions
     363!-------------------------------------------------------------------------------                 
     364      do ind1=1,ngrid
     365
     366      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
     367
     368      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
     369
     370
     371!      zqenv(ind1)=po(ind1)
     372      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
     373      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
     374      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     375      qsatbef=MIN(0.5,qsatbef)
     376      zcor=1./(1.-retv*qsatbef)
     377      qsatbef=qsatbef*zcor
     378      zqsatenv(ind1,ind2)=qsatbef
     379
     380
     381
     382
     383      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
     384      aenv=1./(1.+(alenv*Lv/cppd))
     385      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
     386
     387
     388
     389
     390      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
     391      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
     392      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     393      qsatbef=MIN(0.5,qsatbef)
     394      zcor=1./(1.-retv*qsatbef)
     395      qsatbef=qsatbef*zcor
     396      zqsatth(ind1,ind2)=qsatbef
     397           
     398      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
     399      ath=1./(1.+(alth*Lv/cppd))
     400      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
     401     
     402     
     403
     404!------------------------------------------------------------------------------
     405! Calcul des ?cart-types pour s
     406!------------------------------------------------------------------------------
     407
     408      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
     409      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
     410!       if (paprs(ind1,ind2).gt.90000) then
     411!       ratqs(ind1,ind2)=0.002
     412!       else
     413!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
     414!       endif
     415!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
     416!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
     417!       sigma1s=ratqs(ind1,ind2)*po(ind1)
     418!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
     419 
     420!------------------------------------------------------------------------------
     421! Calcul de l'eau condens?e et de la couverture nuageuse
     422!------------------------------------------------------------------------------
     423      sqrt2pi=sqrt(2.*pi)
     424      xth=sth/(sqrt(2.)*sigma2s)
     425      xenv=senv/(sqrt(2.)*sigma1s)
     426      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
     427      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
     428      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
     429!      ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*cenv(ind1,ind2)
     430
     431
     432
     433      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
     434      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
     435      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
     436!      qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2)
     437     
     438
     439!      print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'
     440
     441
     442!-------------------------------------------------------------------------------
     443!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
     444!-------------------------------------------------------------------------------
     445!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
     446!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
     447      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
     448      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
     449!      deltasenv=aenv*0.01*po(ind1)
     450!     deltasth=ath*0.01*zqta(ind1,ind2)   
     451      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
     452      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
     453      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
     454      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
     455      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
     456      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
     457     
     458      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
     459      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
     460      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
     461
     462      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
     463      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
     464      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
     465      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
     466
     467      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
     468!      qlenv(ind1,ind2)=IntJ
     469!      print*, qlenv(ind1,ind2),'VERIF EAU'
     470
     471
     472      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
     473!      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
     474!      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
     475      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
     476      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
     477      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
     478      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
     479!      qlth(ind1,ind2)=IntJ
     480!      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
     481      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
     482
     483!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     484      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
     485      ctot(ind1,ind2)=0.
     486      qcloud(ind1)=zqsatenv(ind1,ind2)
     487
     488      else
     489               
     490      ctot(ind1,ind2)=ctot(ind1,ind2)
     491      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
     492!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
     493!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
     494
     495      endif 
     496                       
     497     
     498         
     499!     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
     500
     501
     502      else  ! gaussienne environnement seule
     503     
     504      zqenv(ind1)=po(ind1)
     505      Tbef=t(ind1,ind2)
     506      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
     507      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     508      qsatbef=MIN(0.5,qsatbef)
     509      zcor=1./(1.-retv*qsatbef)
     510      qsatbef=qsatbef*zcor
     511      zqsatenv(ind1,ind2)=qsatbef
     512     
     513
     514!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
     515      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
     516      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
     517      aenv=1./(1.+(alenv*Lv/cppd))
     518      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
     519     
     520
     521      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
     522
     523      sqrt2pi=sqrt(2.*pi)
     524      xenv=senv/(sqrt(2.)*sigma1s)
     525      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
     526      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
     527     
     528      if (ctot(ind1,ind2).lt.1.e-3) then
     529      ctot(ind1,ind2)=0.
     530      qcloud(ind1)=zqsatenv(ind1,ind2)
     531
     532      else   
     533               
     534      ctot(ind1,ind2)=ctot(ind1,ind2)
     535      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
     536
     537      endif   
     538 
     539 
     540 
     541 
     542 
     543 
     544      endif   
     545      enddo
     546     
     547      return
     548      end
Note: See TracChangeset for help on using the changeset viewer.