Ignore:
Timestamp:
Mar 5, 2014, 2:19:12 PM (10 years ago)
Author:
lguez
Message:

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

File:
1 moved

Legend:

Unmodified
Added
Removed

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

    r1988 r1992  
    11! $Id$
    2 !
    3       SUBROUTINE nuage (paprs, pplay,
    4      .                  t, pqlwp, pclc, pcltau, pclemi,
    5      .                  pch, pcl, pcm, pct, pctlwp,
    6      e                  ok_aie,
    7      e                  mass_solu_aero, mass_solu_aero_pi,
    8      e                  bl95_b0, bl95_b1,
    9      s                  cldtaupi, re, fl)
    10       USE dimphy
    11       IMPLICIT none
    12 c======================================================================
    13 c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
    14 c Objet: Calculer epaisseur optique et emmissivite des nuages
    15 c======================================================================
    16 c Arguments:
    17 c t-------input-R-temperature
    18 c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
    19 c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
    20 c ok_aie--input-L-apply aerosol indirect effect or not
    21 c mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3]
    22 c mass_solu_aero_pi--input-R-dito, pre-industrial value
    23 c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
    24 c bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
    25 c     
    26 c cldtaupi-output-R-pre-industrial value of cloud optical thickness,
    27 c                   needed for the diagnostics of the aerosol indirect
    28 c                   radiative forcing (see radlwsw)
    29 c re------output-R-Cloud droplet effective radius multiplied by fl [um]
    30 c fl------output-R-Denominator to re, introduced to avoid problems in
    31 c                  the averaging of the output. fl is the fraction of liquid
    32 c                  water clouds within a grid cell     
    33 c
    34 c pcltau--output-R-epaisseur optique des nuages
    35 c pclemi--output-R-emissivite des nuages (0 a 1)
    36 c======================================================================
    37 C
    38 #include "YOMCST.h"
    39 c
    40 cym#include "dimensions.h"
    41 cym#include "dimphy.h"
    42       REAL paprs(klon,klev+1), pplay(klon,klev)
    43       REAL t(klon,klev)
    44 c
    45       REAL pclc(klon,klev)
    46       REAL pqlwp(klon,klev)
    47       REAL pcltau(klon,klev), pclemi(klon,klev)
    48 c
    49       REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
    50 c
    51       LOGICAL lo
    52 c
    53       REAL cetahb, cetamb
    54       PARAMETER (cetahb = 0.45, cetamb = 0.80)
    55 C
    56       INTEGER i, k
    57       REAL zflwp, zradef, zfice, zmsac
    58 c
    59       REAL radius, rad_froid, rad_chaud, rad_chau1, rad_chau2
    60       PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
    61 ccc      PARAMETER (rad_chaud=15.0, rad_froid=35.0)
    62 c sintex initial      PARAMETER (rad_chaud=10.0, rad_froid=30.0)
    63       REAL coef, coef_froi, coef_chau
    64       PARAMETER (coef_chau=0.13, coef_froi=0.09)
    65       REAL seuil_neb, t_glace
    66       PARAMETER (seuil_neb=0.001, t_glace=273.0-15.0)
    67       INTEGER nexpo ! exponentiel pour glace/eau
    68       PARAMETER (nexpo=6)
    69      
    70 cjq for the aerosol indirect effect
    71 cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
    72 cjq     
    73       LOGICAL ok_aie            ! Apply AIE or not?
    74      
    75       REAL mass_solu_aero(klon, klev)    ! total mass concentration for all soluble aerosols[ug m-3]
    76       REAL mass_solu_aero_pi(klon, klev) ! - " - pre-industrial value
    77       REAL cdnc(klon, klev)     ! cloud droplet number concentration [m-3]
    78       REAL re(klon, klev)       ! cloud droplet effective radius [um]
    79       REAL cdnc_pi(klon, klev)     ! cloud droplet number concentration [m-3] (pi value)
    80       REAL re_pi(klon, klev)       ! cloud droplet effective radius [um] (pi value)
    81      
    82       REAL fl(klon, klev)  ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell)
    83      
    84       REAL bl95_b0, bl95_b1     ! Parameter in B&L 95-Formula
    85      
    86       REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
    87 cjq-end     
    88      
    89 ccc      PARAMETER (nexpo=1)
    90 c
    91 c Calculer l'epaisseur optique et l'emmissivite des nuages
    92 c
    93       DO k = 1, klev
    94       DO i = 1, klon
    95          rad_chaud = rad_chau1
    96          IF (k.LE.3) rad_chaud = rad_chau2
    97            
    98          pclc(i,k) = MAX(pclc(i,k), seuil_neb)
    99          zflwp = 1000.*pqlwp(i,k)/RG/pclc(i,k)
    100      .          *(paprs(i,k)-paprs(i,k+1))
    101          zfice = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
    102          zfice = MIN(MAX(zfice,0.0),1.0)
    103          zfice = zfice**nexpo
    104          
    105          IF (ok_aie) THEN
    106             ! Formula "D" of Boucher and Lohmann, Tellus, 1995
    107             !             
    108             cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
    109      .           log(MAX(mass_solu_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3
    110             ! Cloud droplet number concentration (CDNC) is restricted
    111             ! to be within [20, 1000 cm^3]
    112             !
    113             cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
    114             cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
    115      .           log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
    116             cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
    117             !           
    118             !
    119             ! air density: pplay(i,k) / (RD * zT(i,k))
    120             ! factor 1.1: derive effective radius from volume-mean radius
    121             ! factor 1000 is the water density
    122             ! _chaud means that this is the CDR for liquid water clouds
    123             !
    124             rad_chaud =
    125      .           1.1 * ( (pqlwp(i,k) * pplay(i,k) / (RD * T(i,k)) ) 
    126      .               / (4./3. * RPI * 1000. * cdnc(i,k)) )**(1./3.)
    127             !
    128             ! Convert to um. CDR shall be at least 3 um.
    129             !
    130             rad_chaud = MAX(rad_chaud*1.e6, 3.)
    131            
    132             ! For output diagnostics
    133             !
    134             ! Cloud droplet effective radius [um]
    135             !
    136             ! we multiply here with f * xl (fraction of liquid water
    137             ! clouds in the grid cell) to avoid problems in the
    138             ! averaging of the output.
    139             ! In the output of IOIPSL, derive the real cloud droplet
    140             ! effective radius as re/fl
    141             !
    142             fl(i,k) = pclc(i,k)*(1.-zfice)           
    143             re(i,k) = rad_chaud*fl(i,k)
    144            
    145             ! Pre-industrial cloud opt thickness
    146             !
    147             ! "radius" is calculated as rad_chaud above (plus the
    148             ! ice cloud contribution) but using cdnc_pi instead of
    149             ! cdnc.
    150             radius = MAX(1.1e6 * ( (pqlwp(i,k)*pplay(i,k)/(RD*T(i,k))) 
    151      .                / (4./3.*RPI*1000.*cdnc_pi(i,k)) )**(1./3.),
    152      .               3.) * (1.-zfice) + rad_froid * zfice           
    153             cldtaupi(i,k) = 3.0/2.0 * zflwp / radius
    154      .           
    155          ENDIF                  ! ok_aie
    156          
    157          radius = rad_chaud * (1.-zfice) + rad_froid * zfice
    158          coef = coef_chau * (1.-zfice) + coef_froi * zfice
    159          pcltau(i,k) = 3.0/2.0 * zflwp / radius
    160          pclemi(i,k) = 1.0 - EXP( - coef * zflwp)
    161          lo = (pclc(i,k) .LE. seuil_neb)
    162          IF (lo) pclc(i,k) = 0.0
    163          IF (lo) pcltau(i,k) = 0.0
    164          IF (lo) pclemi(i,k) = 0.0
    165          
    166          IF (.NOT.ok_aie) cldtaupi(i,k)=pcltau(i,k)           
    167       ENDDO
    168       ENDDO
    169 ccc      DO k = 1, klev
    170 ccc      DO i = 1, klon
    171 ccc         t(i,k) = t(i,k)
    172 ccc         pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
    173 ccc         lo = pclc(i,k) .GT. (2.*1.e-5)
    174 ccc         zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
    175 ccc     .          /(rg*pclc(i,k))
    176 ccc         zradef = 10.0 + (1.-sigs(k))*45.0
    177 ccc         pcltau(i,k) = 1.5 * zflwp / zradef
    178 ccc         zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
    179 ccc         zmsac = 0.13*(1.0-zfice) + 0.08*zfice
    180 ccc         pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
    181 ccc         if (.NOT.lo) pclc(i,k) = 0.0
    182 ccc         if (.NOT.lo) pcltau(i,k) = 0.0
    183 ccc         if (.NOT.lo) pclemi(i,k) = 0.0
    184 ccc      ENDDO
    185 ccc      ENDDO
    186 cccccc      print*, 'pas de nuage dans le rayonnement'
    187 cccccc      DO k = 1, klev
    188 cccccc      DO i = 1, klon
    189 cccccc         pclc(i,k) = 0.0
    190 cccccc         pcltau(i,k) = 0.0
    191 cccccc         pclemi(i,k) = 0.0
    192 cccccc      ENDDO
    193 cccccc      ENDDO
    194 C
    195 C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
    196 C
    197       DO i = 1, klon
    198          pct(i)=1.0
    199          pch(i)=1.0
    200          pcm(i) = 1.0
    201          pcl(i) = 1.0
    202          pctlwp(i) = 0.0
    203       ENDDO
    204 C
    205       DO k = klev, 1, -1
    206       DO i = 1, klon
    207          pctlwp(i) = pctlwp(i)
    208      .             + pqlwp(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
    209          pct(i) = pct(i)*(1.0-pclc(i,k))
    210          if (pplay(i,k).LE.cetahb*paprs(i,1))
    211      .      pch(i) = pch(i)*(1.0-pclc(i,k))
    212          if (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
    213      .       pplay(i,k).LE.cetamb*paprs(i,1))
    214      .      pcm(i) = pcm(i)*(1.0-pclc(i,k))
    215          if (pplay(i,k).GT.cetamb*paprs(i,1))
    216      .      pcl(i) = pcl(i)*(1.0-pclc(i,k))
    217       ENDDO
    218       ENDDO
    219 C
    220       DO i = 1, klon
    221          pct(i)=1.-pct(i)
    222          pch(i)=1.-pch(i)
    223          pcm(i)=1.-pcm(i)
    224          pcl(i)=1.-pcl(i)
    225       ENDDO
    226 C
    227       RETURN
    228       END
    229       SUBROUTINE diagcld1(paprs,pplay,rain,snow,kbot,ktop,
    230      .                   diafra,dialiq)
    231       use dimphy
    232       IMPLICIT none
    233 c
    234 c Laurent Li (LMD/CNRS), le 12 octobre 1998
    235 c                        (adaptation du code ECMWF)
    236 c
    237 c Dans certains cas, le schema pronostique des nuages n'est
    238 c pas suffisament performant. On a donc besoin de diagnostiquer
    239 c ces nuages. Je dois avouer que c'est une frustration.
    240 c
    241 cym#include "dimensions.h"
    242 cym#include "dimphy.h"
    243 #include "YOMCST.h"
    244 c
    245 c Arguments d'entree:
    246       REAL paprs(klon,klev+1) ! pression (Pa) a inter-couche
    247       REAL pplay(klon,klev) ! pression (Pa) au milieu de couche
    248       REAL t(klon,klev) ! temperature (K)
    249       REAL q(klon,klev) ! humidite specifique (Kg/Kg)
    250       REAL rain(klon) ! pluie convective (kg/m2/s)
    251       REAL snow(klon) ! neige convective (kg/m2/s)
    252       INTEGER ktop(klon) ! sommet de la convection
    253       INTEGER kbot(klon) ! bas de la convection
    254 c
    255 c Arguments de sortie:
    256       REAL diafra(klon,klev) ! fraction nuageuse diagnostiquee
    257       REAL dialiq(klon,klev) ! eau liquide nuageuse
    258 c
    259 c Constantes ajustables:
    260       REAL CANVA, CANVB, CANVH
    261       PARAMETER (CANVA=2.0, CANVB=0.3, CANVH=0.4)
    262       REAL CCA, CCB, CCC
    263       PARAMETER (CCA=0.125, CCB=1.5, CCC=0.8)
    264       REAL CCFCT, CCSCAL
    265       PARAMETER (CCFCT=0.400)
    266       PARAMETER (CCSCAL=1.0E+11)
    267       REAL CETAHB, CETAMB
    268       PARAMETER (CETAHB=0.45, CETAMB=0.80)
    269       REAL CCLWMR
    270       PARAMETER (CCLWMR=1.E-04)
    271       REAL ZEPSCR
    272       PARAMETER (ZEPSCR=1.0E-10)
    273 c
    274 c Variables locales:
    275       INTEGER i, k
    276       REAL zcc(klon)
    277 c
    278 c Initialisation:
    279 c
    280       DO k = 1, klev
    281       DO i = 1, klon
    282          diafra(i,k) = 0.0
    283          dialiq(i,k) = 0.0
    284       ENDDO
    285       ENDDO
    286 c
    287       DO i = 1, klon ! Calculer la fraction nuageuse
    288       zcc(i) = 0.0
    289       IF((rain(i)+snow(i)).GT.0.) THEN
    290          zcc(i)= CCA * LOG(MAX(ZEPSCR,(rain(i)+snow(i))*CCSCAL))-CCB
    291          zcc(i)= MIN(CCC,MAX(0.0,zcc(i)))
    292       ENDIF
    293       ENDDO
    294 c
    295       DO i = 1, klon ! pour traiter les enclumes
    296       diafra(i,ktop(i)) = MAX(diafra(i,ktop(i)),zcc(i)*CCFCT)
    297       IF ((zcc(i).GE.CANVH) .AND.
    298      .    (pplay(i,ktop(i)).LE.CETAHB*paprs(i,1)))
    299      . diafra(i,ktop(i)) = MAX(diafra(i,ktop(i)),
    300      .                         MAX(zcc(i)*CCFCT,CANVA*(zcc(i)-CANVB)))
    301       dialiq(i,ktop(i))=CCLWMR*diafra(i,ktop(i))
    302       ENDDO
    303 c
    304       DO k = 1, klev ! nuages convectifs (sauf enclumes)
    305       DO i = 1, klon
    306       IF (k.LT.ktop(i) .AND. k.GE.kbot(i)) THEN
    307          diafra(i,k)=MAX(diafra(i,k),zcc(i)*CCFCT)
    308          dialiq(i,k)=CCLWMR*diafra(i,k)
    309       ENDIF
    310       ENDDO
    311       ENDDO
    312 c
    313       RETURN
    314       END
    315       SUBROUTINE diagcld2(paprs,pplay,t,q, diafra,dialiq)
    316       use dimphy
    317       IMPLICIT none
    318 c
    319 cym#include "dimensions.h"
    320 cym#include "dimphy.h"
    321 #include "YOMCST.h"
    322 c
    323 c Arguments d'entree:
    324       REAL paprs(klon,klev+1) ! pression (Pa) a inter-couche
    325       REAL pplay(klon,klev) ! pression (Pa) au milieu de couche
    326       REAL t(klon,klev) ! temperature (K)
    327       REAL q(klon,klev) ! humidite specifique (Kg/Kg)
    328 c
    329 c Arguments de sortie:
    330       REAL diafra(klon,klev) ! fraction nuageuse diagnostiquee
    331       REAL dialiq(klon,klev) ! eau liquide nuageuse
    332 c
    333       REAL CETAMB
    334       PARAMETER (CETAMB=0.80)
    335       REAL CLOIA, CLOIB, CLOIC, CLOID
    336       PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.6, CLOID=5.0)
    337 ccc      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
    338       REAL RGAMMAS
    339       PARAMETER (RGAMMAS=0.05)
    340       REAL CRHL
    341       PARAMETER (CRHL=0.15)
    342 ccc      PARAMETER (CRHL=0.70)
    343       REAL t_coup
    344       PARAMETER (t_coup=234.0)
    345 c
    346 c Variables locales:
    347       INTEGER i, k, kb, invb(klon)
    348       REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
    349       REAL zdelta, zcor
    350 c
    351 c Fonctions thermodynamiques:
    352 #include "YOETHF.h"
    353 #include "FCTTRE.h"
    354 c
    355 c Initialisation:
    356 c
    357       DO k = 1, klev
    358       DO i = 1, klon
    359          diafra(i,k) = 0.0
    360          dialiq(i,k) = 0.0
    361       ENDDO
    362       ENDDO
    363 c
    364       DO i = 1, klon
    365          invb(i) = klev
    366          zdthmin(i)=0.0
    367       ENDDO
    368 
    369       DO k = 2, klev/2-1
    370       DO i = 1, klon
    371          zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
    372      .          - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
    373          zdthdp = zdthdp * CLOIA
    374          IF (pplay(i,k).GT.CETAMB*paprs(i,1) .AND.
    375      .       zdthdp.LT.zdthmin(i) ) THEN
    376             zdthmin(i) = zdthdp
    377             invb(i) = k
    378          ENDIF
    379       ENDDO
    380       ENDDO
    381 
    382       DO i = 1, klon
    383          kb=invb(i)
    384          IF (thermcep) THEN
    385            zdelta=MAX(0.,SIGN(1.,RTT-t(i,kb)))
    386            zqs= R2ES*FOEEW(t(i,kb),zdelta)/pplay(i,kb)
    387            zqs=MIN(0.5,zqs)
    388            zcor=1./(1.-RETV*zqs)
    389            zqs=zqs*zcor
    390          ELSE
    391            IF (t(i,kb) .LT. t_coup) THEN
    392               zqs = qsats(t(i,kb)) / pplay(i,kb)
    393            ELSE
    394               zqs = qsatl(t(i,kb)) / pplay(i,kb)
    395            ENDIF
    396          ENDIF
    397          zcll = CLOIB * zdthmin(i) + CLOIC
    398          zcll = MIN(1.0,MAX(0.0,zcll))
    399          zrhb= q(i,kb)/zqs
    400          IF (zcll.GT.0.0.AND.zrhb.LT.CRHL)
    401      .   zcll=zcll*(1.-(CRHL-zrhb)*CLOID)
    402          zcll=MIN(1.0,MAX(0.0,zcll))
    403          diafra(i,kb) = MAX(diafra(i,kb),zcll)
    404          dialiq(i,kb)= diafra(i,kb) * RGAMMAS*zqs
    405       ENDDO
    406 c
    407       RETURN
    408       END
     2
     3SUBROUTINE nuage(paprs, pplay, t, pqlwp, pclc, pcltau, pclemi, pch, pcl, pcm, &
     4    pct, pctlwp, ok_aie, mass_solu_aero, mass_solu_aero_pi, bl95_b0, bl95_b1, &
     5    cldtaupi, re, fl)
     6  USE dimphy
     7  IMPLICIT NONE
     8  ! ======================================================================
     9  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
     10  ! Objet: Calculer epaisseur optique et emmissivite des nuages
     11  ! ======================================================================
     12  ! Arguments:
     13  ! t-------input-R-temperature
     14  ! pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
     15  ! pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
     16  ! ok_aie--input-L-apply aerosol indirect effect or not
     17  ! mass_solu_aero-----input-R-total mass concentration for all soluble
     18  ! aerosols[ug/m^3]
     19  ! mass_solu_aero_pi--input-R-dito, pre-industrial value
     20  ! bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
     21  ! bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
     22
     23  ! cldtaupi-output-R-pre-industrial value of cloud optical thickness,
     24  ! needed for the diagnostics of the aerosol indirect
     25  ! radiative forcing (see radlwsw)
     26  ! re------output-R-Cloud droplet effective radius multiplied by fl [um]
     27  ! fl------output-R-Denominator to re, introduced to avoid problems in
     28  ! the averaging of the output. fl is the fraction of liquid
     29  ! water clouds within a grid cell
     30
     31  ! pcltau--output-R-epaisseur optique des nuages
     32  ! pclemi--output-R-emissivite des nuages (0 a 1)
     33  ! ======================================================================
     34
     35  include "YOMCST.h"
     36
     37  ! ym#include "dimensions.h"
     38  ! ym#include "dimphy.h"
     39  REAL paprs(klon, klev+1), pplay(klon, klev)
     40  REAL t(klon, klev)
     41
     42  REAL pclc(klon, klev)
     43  REAL pqlwp(klon, klev)
     44  REAL pcltau(klon, klev), pclemi(klon, klev)
     45
     46  REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
     47
     48  LOGICAL lo
     49
     50  REAL cetahb, cetamb
     51  PARAMETER (cetahb=0.45, cetamb=0.80)
     52
     53  INTEGER i, k
     54  REAL zflwp, zradef, zfice, zmsac
     55
     56  REAL radius, rad_froid, rad_chaud, rad_chau1, rad_chau2
     57  PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
     58  ! cc      PARAMETER (rad_chaud=15.0, rad_froid=35.0)
     59  ! sintex initial      PARAMETER (rad_chaud=10.0, rad_froid=30.0)
     60  REAL coef, coef_froi, coef_chau
     61  PARAMETER (coef_chau=0.13, coef_froi=0.09)
     62  REAL seuil_neb, t_glace
     63  PARAMETER (seuil_neb=0.001, t_glace=273.0-15.0)
     64  INTEGER nexpo ! exponentiel pour glace/eau
     65  PARAMETER (nexpo=6)
     66
     67  ! jq for the aerosol indirect effect
     68  ! jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
     69  ! jq
     70  LOGICAL ok_aie ! Apply AIE or not?
     71
     72  REAL mass_solu_aero(klon, klev) ! total mass concentration for all soluble aerosols[ug m-3]
     73  REAL mass_solu_aero_pi(klon, klev) ! - " - pre-industrial value
     74  REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
     75  REAL re(klon, klev) ! cloud droplet effective radius [um]
     76  REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)
     77  REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value)
     78
     79  REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds
     80  ! within the grid cell)
     81
     82  REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
     83
     84  REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
     85  ! jq-end
     86
     87  ! cc      PARAMETER (nexpo=1)
     88
     89  ! Calculer l'epaisseur optique et l'emmissivite des nuages
     90
     91  DO k = 1, klev
     92    DO i = 1, klon
     93      rad_chaud = rad_chau1
     94      IF (k<=3) rad_chaud = rad_chau2
     95
     96      pclc(i, k) = max(pclc(i,k), seuil_neb)
     97      zflwp = 1000.*pqlwp(i, k)/rg/pclc(i, k)*(paprs(i,k)-paprs(i,k+1))
     98      zfice = 1.0 - (t(i,k)-t_glace)/(273.13-t_glace)
     99      zfice = min(max(zfice,0.0), 1.0)
     100      zfice = zfice**nexpo
     101
     102      IF (ok_aie) THEN
     103          ! Formula "D" of Boucher and Lohmann, Tellus, 1995
     104          !
     105        cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
     106          1.E-4))/log(10.))*1.E6 !-m-3
     107          ! Cloud droplet number concentration (CDNC) is restricted
     108          ! to be within [20, 1000 cm^3]
     109          !
     110        cdnc(i, k) = min(1000.E6, max(20.E6,cdnc(i,k)))
     111        cdnc_pi(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero_pi(i,k), &
     112          1.E-4))/log(10.))*1.E6 !-m-3
     113        cdnc_pi(i, k) = min(1000.E6, max(20.E6,cdnc_pi(i,k)))
     114          !
     115          !
     116          ! air density: pplay(i,k) / (RD * zT(i,k))
     117          ! factor 1.1: derive effective radius from volume-mean radius
     118          ! factor 1000 is the water density
     119          ! _chaud means that this is the CDR for liquid water clouds
     120          !
     121        rad_chaud = 1.1*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi*1000. &
     122          *cdnc(i,k)))**(1./3.)
     123          !
     124          ! Convert to um. CDR shall be at least 3 um.
     125          !
     126        rad_chaud = max(rad_chaud*1.E6, 3.)
     127
     128          ! For output diagnostics
     129          !
     130          ! Cloud droplet effective radius [um]
     131          !
     132          ! we multiply here with f * xl (fraction of liquid water
     133          ! clouds in the grid cell) to avoid problems in the
     134          ! averaging of the output.
     135          ! In the output of IOIPSL, derive the real cloud droplet
     136          ! effective radius as re/fl
     137          !
     138        fl(i, k) = pclc(i, k)*(1.-zfice)
     139        re(i, k) = rad_chaud*fl(i, k)
     140
     141          ! Pre-industrial cloud opt thickness
     142          !
     143          ! "radius" is calculated as rad_chaud above (plus the
     144          ! ice cloud contribution) but using cdnc_pi instead of
     145          ! cdnc.
     146        radius = max(1.1E6*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi* &
     147          1000.*cdnc_pi(i,k)))**(1./3.), 3.)*(1.-zfice) + rad_froid*zfice
     148        cldtaupi(i, k) = 3.0/2.0*zflwp/radius
     149      END IF ! ok_aie
     150
     151      radius = rad_chaud*(1.-zfice) + rad_froid*zfice
     152      coef = coef_chau*(1.-zfice) + coef_froi*zfice
     153      pcltau(i, k) = 3.0/2.0*zflwp/radius
     154      pclemi(i, k) = 1.0 - exp(-coef*zflwp)
     155      lo = (pclc(i,k)<=seuil_neb)
     156      IF (lo) pclc(i, k) = 0.0
     157      IF (lo) pcltau(i, k) = 0.0
     158      IF (lo) pclemi(i, k) = 0.0
     159
     160      IF (.NOT. ok_aie) cldtaupi(i, k) = pcltau(i, k)
     161    END DO
     162  END DO
     163  ! cc      DO k = 1, klev
     164  ! cc      DO i = 1, klon
     165  ! cc         t(i,k) = t(i,k)
     166  ! cc         pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
     167  ! cc         lo = pclc(i,k) .GT. (2.*1.e-5)
     168  ! cc         zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
     169  ! cc     .          /(rg*pclc(i,k))
     170  ! cc         zradef = 10.0 + (1.-sigs(k))*45.0
     171  ! cc         pcltau(i,k) = 1.5 * zflwp / zradef
     172  ! cc         zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
     173  ! cc         zmsac = 0.13*(1.0-zfice) + 0.08*zfice
     174  ! cc         pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
     175  ! cc         if (.NOT.lo) pclc(i,k) = 0.0
     176  ! cc         if (.NOT.lo) pcltau(i,k) = 0.0
     177  ! cc         if (.NOT.lo) pclemi(i,k) = 0.0
     178  ! cc      ENDDO
     179  ! cc      ENDDO
     180  ! ccccc      print*, 'pas de nuage dans le rayonnement'
     181  ! ccccc      DO k = 1, klev
     182  ! ccccc      DO i = 1, klon
     183  ! ccccc         pclc(i,k) = 0.0
     184  ! ccccc         pcltau(i,k) = 0.0
     185  ! ccccc         pclemi(i,k) = 0.0
     186  ! ccccc      ENDDO
     187  ! ccccc      ENDDO
     188
     189  ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
     190
     191  DO i = 1, klon
     192    pct(i) = 1.0
     193    pch(i) = 1.0
     194    pcm(i) = 1.0
     195    pcl(i) = 1.0
     196    pctlwp(i) = 0.0
     197  END DO
     198
     199  DO k = klev, 1, -1
     200    DO i = 1, klon
     201      pctlwp(i) = pctlwp(i) + pqlwp(i, k)*(paprs(i,k)-paprs(i,k+1))/rg
     202      pct(i) = pct(i)*(1.0-pclc(i,k))
     203      IF (pplay(i,k)<=cetahb*paprs(i,1)) pch(i) = pch(i)*(1.0-pclc(i,k))
     204      IF (pplay(i,k)>cetahb*paprs(i,1) .AND. pplay(i,k)<=cetamb*paprs(i,1)) &
     205        pcm(i) = pcm(i)*(1.0-pclc(i,k))
     206      IF (pplay(i,k)>cetamb*paprs(i,1)) pcl(i) = pcl(i)*(1.0-pclc(i,k))
     207    END DO
     208  END DO
     209
     210  DO i = 1, klon
     211    pct(i) = 1. - pct(i)
     212    pch(i) = 1. - pch(i)
     213    pcm(i) = 1. - pcm(i)
     214    pcl(i) = 1. - pcl(i)
     215  END DO
     216
     217  RETURN
     218END SUBROUTINE nuage
     219SUBROUTINE diagcld1(paprs, pplay, rain, snow, kbot, ktop, diafra, dialiq)
     220  USE dimphy
     221  IMPLICIT NONE
     222
     223  ! Laurent Li (LMD/CNRS), le 12 octobre 1998
     224  ! (adaptation du code ECMWF)
     225
     226  ! Dans certains cas, le schema pronostique des nuages n'est
     227  ! pas suffisament performant. On a donc besoin de diagnostiquer
     228  ! ces nuages. Je dois avouer que c'est une frustration.
     229
     230  ! ym#include "dimensions.h"
     231  ! ym#include "dimphy.h"
     232  include "YOMCST.h"
     233
     234  ! Arguments d'entree:
     235  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
     236  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
     237  REAL t(klon, klev) ! temperature (K)
     238  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
     239  REAL rain(klon) ! pluie convective (kg/m2/s)
     240  REAL snow(klon) ! neige convective (kg/m2/s)
     241  INTEGER ktop(klon) ! sommet de la convection
     242  INTEGER kbot(klon) ! bas de la convection
     243
     244  ! Arguments de sortie:
     245  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
     246  REAL dialiq(klon, klev) ! eau liquide nuageuse
     247
     248  ! Constantes ajustables:
     249  REAL canva, canvb, canvh
     250  PARAMETER (canva=2.0, canvb=0.3, canvh=0.4)
     251  REAL cca, ccb, ccc
     252  PARAMETER (cca=0.125, ccb=1.5, ccc=0.8)
     253  REAL ccfct, ccscal
     254  PARAMETER (ccfct=0.400)
     255  PARAMETER (ccscal=1.0E+11)
     256  REAL cetahb, cetamb
     257  PARAMETER (cetahb=0.45, cetamb=0.80)
     258  REAL cclwmr
     259  PARAMETER (cclwmr=1.E-04)
     260  REAL zepscr
     261  PARAMETER (zepscr=1.0E-10)
     262
     263  ! Variables locales:
     264  INTEGER i, k
     265  REAL zcc(klon)
     266
     267  ! Initialisation:
     268
     269  DO k = 1, klev
     270    DO i = 1, klon
     271      diafra(i, k) = 0.0
     272      dialiq(i, k) = 0.0
     273    END DO
     274  END DO
     275
     276  DO i = 1, klon ! Calculer la fraction nuageuse
     277    zcc(i) = 0.0
     278    IF ((rain(i)+snow(i))>0.) THEN
     279      zcc(i) = cca*log(max(zepscr,(rain(i)+snow(i))*ccscal)) - ccb
     280      zcc(i) = min(ccc, max(0.0,zcc(i)))
     281    END IF
     282  END DO
     283
     284  DO i = 1, klon ! pour traiter les enclumes
     285    diafra(i, ktop(i)) = max(diafra(i,ktop(i)), zcc(i)*ccfct)
     286    IF ((zcc(i)>=canvh) .AND. (pplay(i,ktop(i))<=cetahb*paprs(i, &
     287      1))) diafra(i, ktop(i)) = max(diafra(i,ktop(i)), max(zcc( &
     288      i)*ccfct,canva*(zcc(i)-canvb)))
     289    dialiq(i, ktop(i)) = cclwmr*diafra(i, ktop(i))
     290  END DO
     291
     292  DO k = 1, klev ! nuages convectifs (sauf enclumes)
     293    DO i = 1, klon
     294      IF (k<ktop(i) .AND. k>=kbot(i)) THEN
     295        diafra(i, k) = max(diafra(i,k), zcc(i)*ccfct)
     296        dialiq(i, k) = cclwmr*diafra(i, k)
     297      END IF
     298    END DO
     299  END DO
     300
     301  RETURN
     302END SUBROUTINE diagcld1
     303SUBROUTINE diagcld2(paprs, pplay, t, q, diafra, dialiq)
     304  USE dimphy
     305  IMPLICIT NONE
     306
     307  ! ym#include "dimensions.h"
     308  ! ym#include "dimphy.h"
     309  include "YOMCST.h"
     310
     311  ! Arguments d'entree:
     312  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
     313  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
     314  REAL t(klon, klev) ! temperature (K)
     315  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
     316
     317  ! Arguments de sortie:
     318  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
     319  REAL dialiq(klon, klev) ! eau liquide nuageuse
     320
     321  REAL cetamb
     322  PARAMETER (cetamb=0.80)
     323  REAL cloia, cloib, cloic, cloid
     324  PARAMETER (cloia=1.0E+02, cloib=-10.00, cloic=-0.6, cloid=5.0)
     325  ! cc      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
     326  REAL rgammas
     327  PARAMETER (rgammas=0.05)
     328  REAL crhl
     329  PARAMETER (crhl=0.15)
     330  ! cc      PARAMETER (CRHL=0.70)
     331  REAL t_coup
     332  PARAMETER (t_coup=234.0)
     333
     334  ! Variables locales:
     335  INTEGER i, k, kb, invb(klon)
     336  REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
     337  REAL zdelta, zcor
     338
     339  ! Fonctions thermodynamiques:
     340  include "YOETHF.h"
     341  include "FCTTRE.h"
     342
     343  ! Initialisation:
     344
     345  DO k = 1, klev
     346    DO i = 1, klon
     347      diafra(i, k) = 0.0
     348      dialiq(i, k) = 0.0
     349    END DO
     350  END DO
     351
     352  DO i = 1, klon
     353    invb(i) = klev
     354    zdthmin(i) = 0.0
     355  END DO
     356
     357  DO k = 2, klev/2 - 1
     358    DO i = 1, klon
     359      zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) - &
     360        rd*0.5*(t(i,k)+t(i,k+1))/rcpd/paprs(i, k+1)
     361      zdthdp = zdthdp*cloia
     362      IF (pplay(i,k)>cetamb*paprs(i,1) .AND. zdthdp<zdthmin(i)) THEN
     363        zdthmin(i) = zdthdp
     364        invb(i) = k
     365      END IF
     366    END DO
     367  END DO
     368
     369  DO i = 1, klon
     370    kb = invb(i)
     371    IF (thermcep) THEN
     372      zdelta = max(0., sign(1.,rtt-t(i,kb)))
     373      zqs = r2es*foeew(t(i,kb), zdelta)/pplay(i, kb)
     374      zqs = min(0.5, zqs)
     375      zcor = 1./(1.-retv*zqs)
     376      zqs = zqs*zcor
     377    ELSE
     378      IF (t(i,kb)<t_coup) THEN
     379        zqs = qsats(t(i,kb))/pplay(i, kb)
     380      ELSE
     381        zqs = qsatl(t(i,kb))/pplay(i, kb)
     382      END IF
     383    END IF
     384    zcll = cloib*zdthmin(i) + cloic
     385    zcll = min(1.0, max(0.0,zcll))
     386    zrhb = q(i, kb)/zqs
     387    IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll*(1.-(crhl-zrhb)*cloid)
     388    zcll = min(1.0, max(0.0,zcll))
     389    diafra(i, kb) = max(diafra(i,kb), zcll)
     390    dialiq(i, kb) = diafra(i, kb)*rgammas*zqs
     391  END DO
     392
     393  RETURN
     394END SUBROUTINE diagcld2
Note: See TracChangeset for help on using the changeset viewer.