Ignore:
Timestamp:
Jul 29, 2024, 11:01:04 PM (3 months ago)
Author:
abarral
Message:

Put YOMCST.h into modules

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/branches/Amaury_dev/libf/phylmd/nuage.F90

    r5143 r5144  
    11! $Id$
    22
    3 SUBROUTINE nuage(paprs, pplay, t, pqlwp,picefra, pclc, pcltau, pclemi, pch, pcl, pcm, &
    4     pct, pctlwp, ok_aie, mass_solu_aero, mass_solu_aero_pi, bl95_b0, bl95_b1, distcltop, &
    5     temp_cltop, cldtaupi, re, fl)
     3SUBROUTINE nuage(paprs, pplay, t, pqlwp, picefra, pclc, pcltau, pclemi, pch, pcl, pcm, &
     4        pct, pctlwp, ok_aie, mass_solu_aero, mass_solu_aero_pi, bl95_b0, bl95_b1, distcltop, &
     5        temp_cltop, cldtaupi, re, fl)
    66  USE dimphy
    77  USE lmdz_lscp_tools, ONLY: icefrac_lscp
     
    1111  USE lmdz_clesphys
    1212  USE lmdz_nuage_params ! JBM 3/14
     13  USE lmdz_yomcst
    1314
    1415  IMPLICIT NONE
     
    4142  ! ======================================================================
    4243
    43   include "YOMCST.h"
    44 
    45   REAL paprs(klon, klev+1), pplay(klon, klev)
     44  REAL paprs(klon, klev + 1), pplay(klon, klev)
    4645  REAL t(klon, klev)
    4746
    4847  REAL pclc(klon, klev)
    49   REAL pqlwp(klon, klev), picefra(klon,klev)
     48  REAL pqlwp(klon, klev), picefra(klon, klev)
    5049  REAL pcltau(klon, klev), pclemi(klon, klev)
    5150
    5251  REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
    53   REAL distcltop(klon,klev)
    54   REAL temp_cltop(klon,klev)
     52  REAL distcltop(klon, klev)
     53  REAL temp_cltop(klon, klev)
    5554  LOGICAL lo
    5655
    5756  REAL cetahb, cetamb
    58   PARAMETER (cetahb=0.45, cetamb=0.80)
     57  PARAMETER (cetahb = 0.45, cetamb = 0.80)
    5958
    6059  INTEGER i, k
     
    6261
    6362  REAL radius, rad_chaud
    64 ! JBM (3/14) parameters already defined in nuage.h:
    65 ! REAL rad_froid, rad_chau1, rad_chau2
    66 ! PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
     63  ! JBM (3/14) parameters already defined in nuage.h:
     64  ! REAL rad_froid, rad_chau1, rad_chau2
     65  ! PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
    6766  ! cc      PARAMETER (rad_chaud=15.0, rad_froid=35.0)
    6867  ! sintex initial      PARAMETER (rad_chaud=10.0, rad_froid=30.0)
    6968  REAL coef, coef_froi, coef_chau
    70   PARAMETER (coef_chau=0.13, coef_froi=0.09)
     69  PARAMETER (coef_chau = 0.13, coef_froi = 0.09)
    7170  REAL seuil_neb
    72   PARAMETER (seuil_neb=0.001)
    73 ! JBM (3/14) nexpo is replaced by exposant_glace
    74 ! REAL nexpo ! exponentiel pour glace/eau
    75 ! PARAMETER (nexpo=6.)
     71  PARAMETER (seuil_neb = 0.001)
     72  ! JBM (3/14) nexpo is replaced by exposant_glace
     73  ! REAL nexpo ! exponentiel pour glace/eau
     74  ! PARAMETER (nexpo=6.)
    7675  REAL, PARAMETER :: t_glace_min_old = 258.
    7776  INTEGER, PARAMETER :: exposant_glace_old = 6
     
    104103
    105104  DO k = 1, klev
    106      IF (iflag_t_glace==0) THEN
    107        DO i = 1, klon
    108         zfice(i) = 1.0 - (t(i,k)-t_glace_min_old)/(273.13-t_glace_min_old)
    109         zfice(i) = min(max(zfice(i),0.0), 1.0)
     105    IF (iflag_t_glace==0) THEN
     106      DO i = 1, klon
     107        zfice(i) = 1.0 - (t(i, k) - t_glace_min_old) / (273.13 - t_glace_min_old)
     108        zfice(i) = min(max(zfice(i), 0.0), 1.0)
    110109        zfice(i) = zfice(i)**exposant_glace_old
    111        ENDDO
    112      ELSE ! of IF (iflag_t_glace.EQ.0)
    113 ! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
    114 !       zfice(i) = icefrac_lsc(t(i,k), t_glace_min, &
    115 !                           t_glace_max, exposant_glace)
    116         IF (ok_new_lscp) THEN
    117             CALL icefrac_lscp(klon,t(:,k),iflag_ice_thermo,distcltop(:,k),temp_cltop(:,k),zfice(:),dzfice(:))
    118         ELSE
    119             CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:))
    120 
    121         ENDIF
    122 
    123     IF (ok_new_lscp .AND. ok_icefra_lscp) THEN
    124     ! EV: take the ice fraction directly from the lscp code
    125     ! consistent only for non convective grid points
    126     ! critical for mixed phase clouds
    127         DO i=1,klon
    128         IF (.NOT. ptconv(i,k)) THEN
    129            zfice(i)=picefra(i,k)
    130         ENDIF
     110      ENDDO
     111    ELSE ! of IF (iflag_t_glace.EQ.0)
     112      ! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
     113      !       zfice(i) = icefrac_lsc(t(i,k), t_glace_min, &
     114      !                           t_glace_max, exposant_glace)
     115      IF (ok_new_lscp) THEN
     116        CALL icefrac_lscp(klon, t(:, k), iflag_ice_thermo, distcltop(:, k), temp_cltop(:, k), zfice(:), dzfice(:))
     117      ELSE
     118        CALL icefrac_lsc(klon, t(:, k), pplay(:, k) / paprs(:, 1), zfice(:))
     119
     120      ENDIF
     121
     122      IF (ok_new_lscp .AND. ok_icefra_lscp) THEN
     123        ! EV: take the ice fraction directly from the lscp code
     124        ! consistent only for non convective grid points
     125        ! critical for mixed phase clouds
     126        DO i = 1, klon
     127          IF (.NOT. ptconv(i, k)) THEN
     128            zfice(i) = picefra(i, k)
     129          ENDIF
    131130        ENDDO
    132     ENDIF
    133 
    134 
    135      ENDIF
     131      ENDIF
     132
     133    ENDIF
    136134
    137135    DO i = 1, klon
     
    139137      IF (k<=3) rad_chaud = rad_chau2
    140138
    141       pclc(i, k) = max(pclc(i,k), seuil_neb)
    142       zflwp = 1000.*pqlwp(i, k)/rg/pclc(i, k)*(paprs(i,k)-paprs(i,k+1))
     139      pclc(i, k) = max(pclc(i, k), seuil_neb)
     140      zflwp = 1000. * pqlwp(i, k) / rg / pclc(i, k) * (paprs(i, k) - paprs(i, k + 1))
    143141
    144142      IF (ok_aie) THEN
    145           ! Formula "D" of Boucher and Lohmann, Tellus, 1995
    146 
    147         cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
    148           1.E-4))/log(10.))*1.E6 !-m-3
    149           ! Cloud droplet number concentration (CDNC) is restricted
    150           ! to be within [20, 1000 cm^3]
    151 
    152         cdnc(i, k) = min(1000.E6, max(20.E6,cdnc(i,k)))
    153         cdnc_pi(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero_pi(i,k), &
    154           1.E-4))/log(10.))*1.E6 !-m-3
    155         cdnc_pi(i, k) = min(1000.E6, max(20.E6,cdnc_pi(i,k)))
    156 
    157 
    158           ! air density: pplay(i,k) / (RD * zT(i,k))
    159           ! factor 1.1: derive effective radius from volume-mean radius
    160           ! factor 1000 is the water density
    161           ! _chaud means that this is the CDR for liquid water clouds
    162 
    163         rad_chaud = 1.1*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi*1000. &
    164           *cdnc(i,k)))**(1./3.)
    165 
    166           ! Convert to um. CDR shall be at least 3 um.
    167 
    168         rad_chaud = max(rad_chaud*1.E6, 3.)
    169 
    170           ! For output diagnostics
    171 
    172           ! Cloud droplet effective radius [um]
    173 
    174           ! we multiply here with f * xl (fraction of liquid water
    175           ! clouds in the grid cell) to avoid problems in the
    176           ! averaging of the output.
    177           ! In the output of IOIPSL, derive the real cloud droplet
    178           ! effective radius as re/fl
    179 
    180         fl(i, k) = pclc(i, k)*(1.-zfice(i))
    181         re(i, k) = rad_chaud*fl(i, k)
    182 
    183           ! Pre-industrial cloud opt thickness
    184 
    185           ! "radius" is calculated as rad_chaud above (plus the
    186           ! ice cloud contribution) but using cdnc_pi instead of
    187           ! cdnc.
    188         radius = max(1.1E6*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi* &
    189           1000.*cdnc_pi(i,k)))**(1./3.), 3.)*(1.-zfice(i)) + rad_froid*zfice(i)
    190         cldtaupi(i, k) = 3.0/2.0*zflwp/radius
     143        ! Formula "D" of Boucher and Lohmann, Tellus, 1995
     144
     145        cdnc(i, k) = 10.**(bl95_b0 + bl95_b1 * log(max(mass_solu_aero(i, k), &
     146                1.E-4)) / log(10.)) * 1.E6 !-m-3
     147        ! Cloud droplet number concentration (CDNC) is restricted
     148        ! to be within [20, 1000 cm^3]
     149
     150        cdnc(i, k) = min(1000.E6, max(20.E6, cdnc(i, k)))
     151        cdnc_pi(i, k) = 10.**(bl95_b0 + bl95_b1 * log(max(mass_solu_aero_pi(i, k), &
     152                1.E-4)) / log(10.)) * 1.E6 !-m-3
     153        cdnc_pi(i, k) = min(1000.E6, max(20.E6, cdnc_pi(i, k)))
     154
     155
     156        ! air density: pplay(i,k) / (RD * zT(i,k))
     157        ! factor 1.1: derive effective radius from volume-mean radius
     158        ! factor 1000 is the water density
     159        ! _chaud means that this is the CDR for liquid water clouds
     160
     161        rad_chaud = 1.1 * ((pqlwp(i, k) * pplay(i, k) / (rd * t(i, k))) / (4. / 3. * rpi * 1000. &
     162                * cdnc(i, k)))**(1. / 3.)
     163
     164        ! Convert to um. CDR shall be at least 3 um.
     165
     166        rad_chaud = max(rad_chaud * 1.E6, 3.)
     167
     168        ! For output diagnostics
     169
     170        ! Cloud droplet effective radius [um]
     171
     172        ! we multiply here with f * xl (fraction of liquid water
     173        ! clouds in the grid cell) to avoid problems in the
     174        ! averaging of the output.
     175        ! In the output of IOIPSL, derive the real cloud droplet
     176        ! effective radius as re/fl
     177
     178        fl(i, k) = pclc(i, k) * (1. - zfice(i))
     179        re(i, k) = rad_chaud * fl(i, k)
     180
     181        ! Pre-industrial cloud opt thickness
     182
     183        ! "radius" is calculated as rad_chaud above (plus the
     184        ! ice cloud contribution) but using cdnc_pi instead of
     185        ! cdnc.
     186        radius = max(1.1E6 * ((pqlwp(i, k) * pplay(i, k) / (rd * t(i, k))) / (4. / 3. * rpi * &
     187                1000. * cdnc_pi(i, k)))**(1. / 3.), 3.) * (1. - zfice(i)) + rad_froid * zfice(i)
     188        cldtaupi(i, k) = 3.0 / 2.0 * zflwp / radius
    191189      END IF ! ok_aie
    192190
    193       radius = rad_chaud*(1.-zfice(i)) + rad_froid*zfice(i)
    194       coef = coef_chau*(1.-zfice(i)) + coef_froi*zfice(i)
    195       pcltau(i, k) = 3.0/2.0*zflwp/radius
    196       pclemi(i, k) = 1.0 - exp(-coef*zflwp)
    197       lo = (pclc(i,k)<=seuil_neb)
     191      radius = rad_chaud * (1. - zfice(i)) + rad_froid * zfice(i)
     192      coef = coef_chau * (1. - zfice(i)) + coef_froi * zfice(i)
     193      pcltau(i, k) = 3.0 / 2.0 * zflwp / radius
     194      pclemi(i, k) = 1.0 - exp(-coef * zflwp)
     195      lo = (pclc(i, k)<=seuil_neb)
    198196      IF (lo) pclc(i, k) = 0.0
    199197      IF (lo) pcltau(i, k) = 0.0
     
    241239  DO k = klev, 1, -1
    242240    DO i = 1, klon
    243       pctlwp(i) = pctlwp(i) + pqlwp(i, k)*(paprs(i,k)-paprs(i,k+1))/rg
    244       pct(i) = pct(i)*(1.0-pclc(i,k))
    245       IF (pplay(i,k)<=cetahb*paprs(i,1)) pch(i) = pch(i)*(1.0-pclc(i,k))
    246       IF (pplay(i,k)>cetahb*paprs(i,1) .AND. pplay(i,k)<=cetamb*paprs(i,1)) &
    247         pcm(i) = pcm(i)*(1.0-pclc(i,k))
    248       IF (pplay(i,k)>cetamb*paprs(i,1)) pcl(i) = pcl(i)*(1.0-pclc(i,k))
     241      pctlwp(i) = pctlwp(i) + pqlwp(i, k) * (paprs(i, k) - paprs(i, k + 1)) / rg
     242      pct(i) = pct(i) * (1.0 - pclc(i, k))
     243      IF (pplay(i, k)<=cetahb * paprs(i, 1)) pch(i) = pch(i) * (1.0 - pclc(i, k))
     244      IF (pplay(i, k)>cetahb * paprs(i, 1) .AND. pplay(i, k)<=cetamb * paprs(i, 1)) &
     245              pcm(i) = pcm(i) * (1.0 - pclc(i, k))
     246      IF (pplay(i, k)>cetamb * paprs(i, 1)) pcl(i) = pcl(i) * (1.0 - pclc(i, k))
    249247    END DO
    250248  END DO
     
    257255  END DO
    258256
    259 
    260257END SUBROUTINE nuage
    261258SUBROUTINE diagcld1(paprs, pplay, rain, snow, kbot, ktop, diafra, dialiq)
    262259  USE dimphy
     260  USE lmdz_yomcst
     261
    263262  IMPLICIT NONE
    264263
     
    270269  ! ces nuages. Je dois avouer que c'est une frustration.
    271270
    272   include "YOMCST.h"
    273 
    274271  ! Arguments d'entree:
    275   REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
     272  REAL paprs(klon, klev + 1) ! pression (Pa) a inter-couche
    276273  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
    277274  REAL t(klon, klev) ! temperature (K)
     
    288285  ! Constantes ajustables:
    289286  REAL canva, canvb, canvh
    290   PARAMETER (canva=2.0, canvb=0.3, canvh=0.4)
     287  PARAMETER (canva = 2.0, canvb = 0.3, canvh = 0.4)
    291288  REAL cca, ccb, ccc
    292   PARAMETER (cca=0.125, ccb=1.5, ccc=0.8)
     289  PARAMETER (cca = 0.125, ccb = 1.5, ccc = 0.8)
    293290  REAL ccfct, ccscal
    294   PARAMETER (ccfct=0.400)
    295   PARAMETER (ccscal=1.0E+11)
     291  PARAMETER (ccfct = 0.400)
     292  PARAMETER (ccscal = 1.0E+11)
    296293  REAL cetahb, cetamb
    297   PARAMETER (cetahb=0.45, cetamb=0.80)
     294  PARAMETER (cetahb = 0.45, cetamb = 0.80)
    298295  REAL cclwmr
    299   PARAMETER (cclwmr=1.E-04)
     296  PARAMETER (cclwmr = 1.E-04)
    300297  REAL zepscr
    301   PARAMETER (zepscr=1.0E-10)
     298  PARAMETER (zepscr = 1.0E-10)
    302299
    303300  ! Variables locales:
     
    316313  DO i = 1, klon ! Calculer la fraction nuageuse
    317314    zcc(i) = 0.0
    318     IF ((rain(i)+snow(i))>0.) THEN
    319       zcc(i) = cca*log(max(zepscr,(rain(i)+snow(i))*ccscal)) - ccb
    320       zcc(i) = min(ccc, max(0.0,zcc(i)))
     315    IF ((rain(i) + snow(i))>0.) THEN
     316      zcc(i) = cca * log(max(zepscr, (rain(i) + snow(i)) * ccscal)) - ccb
     317      zcc(i) = min(ccc, max(0.0, zcc(i)))
    321318    END IF
    322319  END DO
    323320
    324321  DO i = 1, klon ! pour traiter les enclumes
    325     diafra(i, ktop(i)) = max(diafra(i,ktop(i)), zcc(i)*ccfct)
    326     IF ((zcc(i)>=canvh) .AND. (pplay(i,ktop(i))<=cetahb*paprs(i, &
    327       1))) diafra(i, ktop(i)) = max(diafra(i,ktop(i)), max(zcc( &
    328       i)*ccfct,canva*(zcc(i)-canvb)))
    329     dialiq(i, ktop(i)) = cclwmr*diafra(i, ktop(i))
     322    diafra(i, ktop(i)) = max(diafra(i, ktop(i)), zcc(i) * ccfct)
     323    IF ((zcc(i)>=canvh) .AND. (pplay(i, ktop(i))<=cetahb * paprs(i, &
     324            1))) diafra(i, ktop(i)) = max(diafra(i, ktop(i)), max(zcc(&
     325            i) * ccfct, canva * (zcc(i) - canvb)))
     326    dialiq(i, ktop(i)) = cclwmr * diafra(i, ktop(i))
    330327  END DO
    331328
     
    333330    DO i = 1, klon
    334331      IF (k<ktop(i) .AND. k>=kbot(i)) THEN
    335         diafra(i, k) = max(diafra(i,k), zcc(i)*ccfct)
    336         dialiq(i, k) = cclwmr*diafra(i, k)
     332        diafra(i, k) = max(diafra(i, k), zcc(i) * ccfct)
     333        dialiq(i, k) = cclwmr * diafra(i, k)
    337334      END IF
    338335    END DO
    339336  END DO
    340 
    341337
    342338END SUBROUTINE diagcld1
    343339SUBROUTINE diagcld2(paprs, pplay, t, q, diafra, dialiq)
    344340  USE dimphy
    345   USE lmdz_YOETHF
     341  USE lmdz_yoethf
    346342  USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
     343  USE lmdz_yomcst
    347344
    348345  IMPLICIT NONE
    349346
    350   include "YOMCST.h"
    351 
    352347  ! Arguments d'entree:
    353   REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
     348  REAL paprs(klon, klev + 1) ! pression (Pa) a inter-couche
    354349  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
    355350  REAL t(klon, klev) ! temperature (K)
     
    361356
    362357  REAL cetamb
    363   PARAMETER (cetamb=0.80)
     358  PARAMETER (cetamb = 0.80)
    364359  REAL cloia, cloib, cloic, cloid
    365   PARAMETER (cloia=1.0E+02, cloib=-10.00, cloic=-0.6, cloid=5.0)
     360  PARAMETER (cloia = 1.0E+02, cloib = -10.00, cloic = -0.6, cloid = 5.0)
    366361  ! cc      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
    367362  REAL rgammas
    368   PARAMETER (rgammas=0.05)
     363  PARAMETER (rgammas = 0.05)
    369364  REAL crhl
    370   PARAMETER (crhl=0.15)
     365  PARAMETER (crhl = 0.15)
    371366  ! cc      PARAMETER (CRHL=0.70)
    372367  REAL t_coup
    373   PARAMETER (t_coup=234.0)
     368  PARAMETER (t_coup = 234.0)
    374369
    375370  ! Variables locales:
     
    392387  END DO
    393388
    394   DO k = 2, klev/2 - 1
    395     DO i = 1, klon
    396       zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) - &
    397         rd*0.5*(t(i,k)+t(i,k+1))/rcpd/paprs(i, k+1)
    398       zdthdp = zdthdp*cloia
    399       IF (pplay(i,k)>cetamb*paprs(i,1) .AND. zdthdp<zdthmin(i)) THEN
     389  DO k = 2, klev / 2 - 1
     390    DO i = 1, klon
     391      zdthdp = (t(i, k) - t(i, k + 1)) / (pplay(i, k) - pplay(i, k + 1)) - &
     392              rd * 0.5 * (t(i, k) + t(i, k + 1)) / rcpd / paprs(i, k + 1)
     393      zdthdp = zdthdp * cloia
     394      IF (pplay(i, k)>cetamb * paprs(i, 1) .AND. zdthdp<zdthmin(i)) THEN
    400395        zdthmin(i) = zdthdp
    401396        invb(i) = k
     
    407402    kb = invb(i)
    408403    IF (thermcep) THEN
    409       zdelta = max(0., sign(1.,rtt-t(i,kb)))
    410       zqs = r2es*foeew(t(i,kb), zdelta)/pplay(i, kb)
     404      zdelta = max(0., sign(1., rtt - t(i, kb)))
     405      zqs = r2es * foeew(t(i, kb), zdelta) / pplay(i, kb)
    411406      zqs = min(0.5, zqs)
    412       zcor = 1./(1.-retv*zqs)
    413       zqs = zqs*zcor
     407      zcor = 1. / (1. - retv * zqs)
     408      zqs = zqs * zcor
    414409    ELSE
    415       IF (t(i,kb)<t_coup) THEN
    416         zqs = qsats(t(i,kb))/pplay(i, kb)
     410      IF (t(i, kb)<t_coup) THEN
     411        zqs = qsats(t(i, kb)) / pplay(i, kb)
    417412      ELSE
    418         zqs = qsatl(t(i,kb))/pplay(i, kb)
     413        zqs = qsatl(t(i, kb)) / pplay(i, kb)
    419414      END IF
    420415    END IF
    421     zcll = cloib*zdthmin(i) + cloic
    422     zcll = min(1.0, max(0.0,zcll))
    423     zrhb = q(i, kb)/zqs
    424     IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll*(1.-(crhl-zrhb)*cloid)
    425     zcll = min(1.0, max(0.0,zcll))
    426     diafra(i, kb) = max(diafra(i,kb), zcll)
    427     dialiq(i, kb) = diafra(i, kb)*rgammas*zqs
    428   END DO
    429 
     416    zcll = cloib * zdthmin(i) + cloic
     417    zcll = min(1.0, max(0.0, zcll))
     418    zrhb = q(i, kb) / zqs
     419    IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll * (1. - (crhl - zrhb) * cloid)
     420    zcll = min(1.0, max(0.0, zcll))
     421    diafra(i, kb) = max(diafra(i, kb), zcll)
     422    dialiq(i, kb) = diafra(i, kb) * rgammas * zqs
     423  END DO
    430424
    431425END SUBROUTINE diagcld2
Note: See TracChangeset for help on using the changeset viewer.