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/fisrtilp_tr.F90

    r5143 r5144  
    1 
    21! $Id$
    32
    43
    54SUBROUTINE fisrtilp_tr(dtime, paprs, pplay, t, q, ratqs, d_t, d_q, d_ql, &
    6     rneb, radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, &
    7     frac_nucl, prfl, psfl, rhcl) ! relative humidity in clear sky (needed for aer optical
     5        rneb, radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, &
     6        frac_nucl, prfl, psfl, rhcl) ! relative humidity in clear sky (needed for aer optical
    87  ! properties; aeropt.F)
    9 
    108
    119  USE dimphy
    1210  USE lmdz_print_control, ONLY: lunout
    13   USE lmdz_YOETHF
     11  USE lmdz_yoethf
    1412  USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
     13  USE lmdz_yomcst
    1514
    1615  IMPLICIT NONE
    … …  
    2221  ! ======================================================================
    2322  ! ======================================================================
    24   include "YOMCST.h"
    2523
    2624  ! Arguments:
    2725
    2826  REAL dtime ! intervalle du temps (s)
    29   REAL paprs(klon, klev+1) ! pression a inter-couche
     27  REAL paprs(klon, klev + 1) ! pression a inter-couche
    3028  REAL pplay(klon, klev) ! pression au milieu de couche
    3129  REAL t(klon, klev) ! temperature (K)
    … …  
    3836  REAL rain(klon) ! pluies (mm/s)
    3937  REAL snow(klon) ! neige (mm/s)
    40   REAL prfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
    41   REAL psfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
     38  REAL prfl(klon, klev + 1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
     39  REAL psfl(klon, klev + 1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
    4240
    4341  ! jq   For aerosol opt properties needed (see aeropt.F)
    … …  
    6159
    6260  REAL seuil_neb ! un nuage existe vraiment au-dela
    63   PARAMETER (seuil_neb=0.001)
     61  PARAMETER (seuil_neb = 0.001)
    6462  REAL ct ! inverse du temps pour qu'un nuage precipite
    65   PARAMETER (ct=1./1800.)
     63  PARAMETER (ct = 1. / 1800.)
    6664  REAL cl ! seuil de precipitation
    67   PARAMETER (cl=2.6E-4)
     65  PARAMETER (cl = 2.6E-4)
    6866  ! cc      PARAMETER (cl=2.3e-4)
    6967  ! cc      PARAMETER (cl=2.0e-4)
    7068  INTEGER ninter ! sous-intervals pour la precipitation
    71   PARAMETER (ninter=5)
     69  PARAMETER (ninter = 5)
    7270  LOGICAL evap_prec ! evaporation de la pluie
    73   PARAMETER (evap_prec=.TRUE.)
     71  PARAMETER (evap_prec = .TRUE.)
    7472  REAL coef_eva
    75   PARAMETER (coef_eva=2.0E-05)
     73  PARAMETER (coef_eva = 2.0E-05)
    7674  LOGICAL calcrat ! calculer ratqs au lieu de fixer sa valeur
    7775  REAL ratqs(klon, klev) ! determine la largeur de distribution de vapeur
    78   PARAMETER (calcrat=.TRUE.)
     76  PARAMETER (calcrat = .TRUE.)
    7977  REAL zx_min, rat_max
    80   PARAMETER (zx_min=1.0, rat_max=0.01)
     78  PARAMETER (zx_min = 1.0, rat_max = 0.01)
    8179  REAL zx_max, rat_min
    82   PARAMETER (zx_max=0.1, rat_min=0.3)
     80  PARAMETER (zx_max = 0.1, rat_min = 0.3)
    8381  REAL zx
    8482
    8583  LOGICAL cpartiel ! condensation partielle
    86   PARAMETER (cpartiel=.TRUE.)
     84  PARAMETER (cpartiel = .TRUE.)
    8785  REAL t_coup
    88   PARAMETER (t_coup=234.0)
     86  PARAMETER (t_coup = 234.0)
    8987
    9088  ! Variables locales:
    … …  
    125123  REAL fallv ! vitesse de chute pour crystaux de glace
    126124  REAL zzz
    127   fallv(zzz) = 3.29/2.0*((zzz)**0.16)
     125  DATA appel1er/.TRUE./
     126  fallv(zzz) = 3.29 / 2.0 * ((zzz)**0.16)
    128127  ! cc      fallv (zzz) = 3.29/3.0 * ((zzz)**0.16)
    129128  ! cc      fallv (zzz) = 3.29 * ((zzz)**0.16)
    130 
    131   DATA appel1er/.TRUE./
    132129
    133130  IF (appel1er) THEN
    … …  
    137134    WRITE (lunout, *) 'fisrtilp, evap_prec:', evap_prec
    138135    WRITE (lunout, *) 'fisrtilp, cpartiel:', cpartiel
    139     IF (abs(dtime/real(ninter)-360.0)>0.001) THEN
     136    IF (abs(dtime / real(ninter) - 360.0)>0.001) THEN
    140137      WRITE (lunout, *) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
    141138      WRITE (lunout, *) 'Je prefere un sous-intervalle de 6 minutes'
    … …  
    229226        IF (zrfl(i)>0.) THEN
    230227          IF (thermcep) THEN
    231             zdelta = max(0., sign(1.,rtt-zt(i)))
    232             zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
     228            zdelta = max(0., sign(1., rtt - zt(i)))
     229            zqs(i) = r2es * foeew(zt(i), zdelta) / pplay(i, k)
    233230            zqs(i) = min(0.5, zqs(i))
    234             zcor = 1./(1.-retv*zqs(i))
    235             zqs(i) = zqs(i)*zcor
     231            zcor = 1. / (1. - retv * zqs(i))
     232            zqs(i) = zqs(i) * zcor
    236233          ELSE
    237234            IF (zt(i)<t_coup) THEN
    238               zqs(i) = qsats(zt(i))/pplay(i, k)
     235              zqs(i) = qsats(zt(i)) / pplay(i, k)
    239236            ELSE
    240               zqs(i) = qsatl(zt(i))/pplay(i, k)
     237              zqs(i) = qsatl(zt(i)) / pplay(i, k)
    241238            END IF
    242239          END IF
    243           zqev = max(0.0, (zqs(i)-zq(i))*zneb(i))
    244           zqevt = coef_eva*(1.0-zq(i)/zqs(i))*sqrt(zrfl(i))* &
    245             (paprs(i,k)-paprs(i,k+1))/pplay(i, k)*zt(i)*rd/rg
    246           zqevt = max(0.0, min(zqevt,zrfl(i)))*rg*dtime/ &
    247             (paprs(i,k)-paprs(i,k+1))
     240          zqev = max(0.0, (zqs(i) - zq(i)) * zneb(i))
     241          zqevt = coef_eva * (1.0 - zq(i) / zqs(i)) * sqrt(zrfl(i)) * &
     242                  (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * rd / rg
     243          zqevt = max(0.0, min(zqevt, zrfl(i))) * rg * dtime / &
     244                  (paprs(i, k) - paprs(i, k + 1))
    248245          zqev = min(zqev, zqevt)
    249           zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))/rg/dtime
    250           zq(i) = zq(i) - (zrfln(i)-zrfl(i))*(rg/(paprs(i,k)-paprs(i, &
    251             k+1)))*dtime
    252           zt(i) = zt(i) + (zrfln(i)-zrfl(i))*(rg/(paprs(i,k)-paprs(i, &
    253             k+1)))*dtime*rlvtt/rcpd/(1.0+rvtmp2*zq(i))
     246          zrfln(i) = zrfl(i) - zqev * (paprs(i, k) - paprs(i, k + 1)) / rg / dtime
     247          zq(i) = zq(i) - (zrfln(i) - zrfl(i)) * (rg / (paprs(i, k) - paprs(i, &
     248                  k + 1))) * dtime
     249          zt(i) = zt(i) + (zrfln(i) - zrfl(i)) * (rg / (paprs(i, k) - paprs(i, &
     250                  k + 1))) * dtime * rlvtt / rcpd / (1.0 + rvtmp2 * zq(i))
    254251          zrfl(i) = zrfln(i)
    255252        END IF
    … …  
    261258    IF (thermcep) THEN
    262259      DO i = 1, klon
    263         zdelta = max(0., sign(1.,rtt-zt(i)))
    264         zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta
    265         zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i))
    266         zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
     260        zdelta = max(0., sign(1., rtt - zt(i)))
     261        zcvm5 = r5les * rlvtt * (1. - zdelta) + r5ies * rlstt * zdelta
     262        zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * zq(i))
     263        zqs(i) = r2es * foeew(zt(i), zdelta) / pplay(i, k)
    267264        zqs(i) = min(0.5, zqs(i))
    268         zcor = 1./(1.-retv*zqs(i))
    269         zqs(i) = zqs(i)*zcor
     265        zcor = 1. / (1. - retv * zqs(i))
     266        zqs(i) = zqs(i) * zcor
    270267        zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)
    271268      END DO
    … …  
    273270      DO i = 1, klon
    274271        IF (zt(i)<t_coup) THEN
    275           zqs(i) = qsats(zt(i))/pplay(i, k)
     272          zqs(i) = qsats(zt(i)) / pplay(i, k)
    276273          zdqs(i) = dqsats(zt(i), zqs(i))
    277274        ELSE
    278           zqs(i) = qsatl(zt(i))/pplay(i, k)
     275          zqs(i) = qsatl(zt(i)) / pplay(i, k)
    279276          zdqs(i) = dqsatl(zt(i), zqs(i))
    280277        END IF
    … …  
    288285      DO i = 1, klon
    289286
    290         zdelq = ratqs(i, k)*zq(i)
    291         rneb(i, k) = (zq(i)+zdelq-zqs(i))/(2.0*zdelq)
    292         zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
    293         IF (rneb(i,k)<=0.0) zqn(i) = 0.0
    294         IF (rneb(i,k)>=1.0) zqn(i) = zq(i)
    295         rneb(i, k) = max(0.0, min(1.0,rneb(i,k)))
    296         zcond(i) = max(0.0, zqn(i)-zqs(i))*rneb(i, k)/(1.+zdqs(i))
     287        zdelq = ratqs(i, k) * zq(i)
     288        rneb(i, k) = (zq(i) + zdelq - zqs(i)) / (2.0 * zdelq)
     289        zqn(i) = (zq(i) + zdelq + zqs(i)) / 2.0
     290        IF (rneb(i, k)<=0.0) zqn(i) = 0.0
     291        IF (rneb(i, k)>=1.0) zqn(i) = zq(i)
     292        rneb(i, k) = max(0.0, min(1.0, rneb(i, k)))
     293        zcond(i) = max(0.0, zqn(i) - zqs(i)) * rneb(i, k) / (1. + zdqs(i))
    297294
    298295        ! --Olivier
    299         rhcl(i, k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i)
    300         IF (rneb(i,k)<=0.0) rhcl(i, k) = zq(i)/zqs(i)
    301         IF (rneb(i,k)>=1.0) rhcl(i, k) = 1.0
     296        rhcl(i, k) = (zqs(i) + zq(i) - zdelq) / 2. / zqs(i)
     297        IF (rneb(i, k)<=0.0) rhcl(i, k) = zq(i) / zqs(i)
     298        IF (rneb(i, k)>=1.0) rhcl(i, k) = 1.0
    302299        ! --fin
    303300
    … …  
    310307          rneb(i, k) = 0.0
    311308        END IF
    312         zcond(i) = max(0.0, zq(i)-zqs(i))/(1.+zdqs(i))
     309        zcond(i) = max(0.0, zq(i) - zqs(i)) / (1. + zdqs(i))
    313310      END DO
    314311    END IF
    … …  
    316313    DO i = 1, klon
    317314      zq(i) = zq(i) - zcond(i)
    318       zt(i) = zt(i) + zcond(i)*rlvtt/rcpd
     315      zt(i) = zt(i) + zcond(i) * rlvtt / rcpd
    319316    END DO
    320317
    … …  
    322319
    323320    DO i = 1, klon
    324       IF (rneb(i,k)>0.0) THEN
     321      IF (rneb(i, k)>0.0) THEN
    325322        zoliq(i) = zcond(i)
    326         zrho(i) = pplay(i, k)/zt(i)/rd
    327         zdz(i) = (paprs(i,k)-paprs(i,k+1))/(zrho(i)*rg)
    328         zfice(i) = 1.0 - (zt(i)-ztglace)/(273.13-ztglace)
    329         zfice(i) = min(max(zfice(i),0.0), 1.0)
     323        zrho(i) = pplay(i, k) / zt(i) / rd
     324        zdz(i) = (paprs(i, k) - paprs(i, k + 1)) / (zrho(i) * rg)
     325        zfice(i) = 1.0 - (zt(i) - ztglace) / (273.13 - ztglace)
     326        zfice(i) = min(max(zfice(i), 0.0), 1.0)
    330327        zfice(i) = zfice(i)**nexpo
    331         zneb(i) = max(rneb(i,k), seuil_neb)
    332         radliq(i, k) = zoliq(i)/real(ninter+1)
     328        zneb(i) = max(rneb(i, k), seuil_neb)
     329        radliq(i, k) = zoliq(i) / real(ninter + 1)
    333330      END IF
    334331    END DO
    … …  
    336333    DO n = 1, ninter
    337334      DO i = 1, klon
    338         IF (rneb(i,k)>0.0) THEN
    339           zchau(i) = ct*dtime/real(ninter)*zoliq(i)* &
    340             (1.0-exp(-(zoliq(i)/zneb(i)/cl)**2))*(1.-zfice(i))
    341           zrhol(i) = zrho(i)*zoliq(i)/zneb(i)
    342           zfroi(i) = dtime/real(ninter)/zdz(i)*zoliq(i)*fallv(zrhol(i))* &
    343             zfice(i)
     335        IF (rneb(i, k)>0.0) THEN
     336          zchau(i) = ct * dtime / real(ninter) * zoliq(i) * &
     337                  (1.0 - exp(-(zoliq(i) / zneb(i) / cl)**2)) * (1. - zfice(i))
     338          zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
     339          zfroi(i) = dtime / real(ninter) / zdz(i) * zoliq(i) * fallv(zrhol(i)) * &
     340                  zfice(i)
    344341          ztot(i) = zchau(i) + zfroi(i)
    345342          IF (zneb(i)==seuil_neb) ztot(i) = 0.0
    346           ztot(i) = min(max(ztot(i),0.0), zoliq(i))
    347           zoliq(i) = max(zoliq(i)-ztot(i), 0.0)
    348           radliq(i, k) = radliq(i, k) + zoliq(i)/real(ninter+1)
    349         END IF
    350       END DO
    351     END DO
    352 
    353     DO i = 1, klon
    354       IF (rneb(i,k)>0.0) THEN
     343          ztot(i) = min(max(ztot(i), 0.0), zoliq(i))
     344          zoliq(i) = max(zoliq(i) - ztot(i), 0.0)
     345          radliq(i, k) = radliq(i, k) + zoliq(i) / real(ninter + 1)
     346        END IF
     347      END DO
     348    END DO
     349
     350    DO i = 1, klon
     351      IF (rneb(i, k)>0.0) THEN
    355352        d_ql(i, k) = zoliq(i)
    356         zrfl(i) = zrfl(i) + max(zcond(i)-zoliq(i), 0.0)*(paprs(i,k)-paprs(i,k &
    357           +1))/(rg*dtime)
     353        zrfl(i) = zrfl(i) + max(zcond(i) - zoliq(i), 0.0) * (paprs(i, k) - paprs(i, k &
     354                + 1)) / (rg * dtime)
    358355      END IF
    359356      IF (zt(i)<rtt) THEN
    … …  
    375372    DO i = 1, klon
    376373
    377       zprec_cond(i) = max(zcond(i)-zoliq(i), 0.0)*(paprs(i,k)-paprs(i,k+1))/ &
    378         rg
    379       IF (rneb(i,k)>0.0 .AND. zprec_cond(i)>0.) THEN
     374      zprec_cond(i) = max(zcond(i) - zoliq(i), 0.0) * (paprs(i, k) - paprs(i, k + 1)) / &
     375              rg
     376      IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN
    380377        ! AA lessivage nucleation LMD5 dans la couche elle-meme
    381         IF (t(i,k)>=ztglace) THEN
     378        IF (t(i, k)>=ztglace) THEN
    382379          zalpha_tr = a_tr_sca(3)
    383380        ELSE
    384381          zalpha_tr = a_tr_sca(4)
    385382        END IF
    386         zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
    387         pfrac_nucl(i, k) = pfrac_nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
    388         frac_nucl(i, k) = 1. - zneb(i)*zfrac_lessi
     383        zfrac_lessi = 1. - exp(zalpha_tr * zprec_cond(i) / zneb(i))
     384        pfrac_nucl(i, k) = pfrac_nucl(i, k) * (1. - zneb(i) * zfrac_lessi)
     385        frac_nucl(i, k) = 1. - zneb(i) * zfrac_lessi
    389386
    390387        ! nucleation avec un facteur -1 au lieu de -0.5
    391         zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))
    392         pfrac_1nucl(i, k) = pfrac_1nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
     388        zfrac_lessi = 1. - exp(-zprec_cond(i) / zneb(i))
     389        pfrac_1nucl(i, k) = pfrac_1nucl(i, k) * (1. - zneb(i) * zfrac_lessi)
    393390      END IF
    394391
    … …  
    398395    DO kk = k - 1, 1, -1
    399396      DO i = 1, klon
    400         IF (rneb(i,k)>0.0 .AND. zprec_cond(i)>0.) THEN
    401           IF (t(i,kk)>=ztglace) THEN
     397        IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN
     398          IF (t(i, kk)>=ztglace) THEN
    402399            zalpha_tr = a_tr_sca(1)
    403400          ELSE
    404401            zalpha_tr = a_tr_sca(2)
    405402          END IF
    406           zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
    407           pfrac_impa(i, kk) = pfrac_impa(i, kk)*(1.-zneb(i)*zfrac_lessi)
    408           frac_impa(i, kk) = 1. - zneb(i)*zfrac_lessi
     403          zfrac_lessi = 1. - exp(zalpha_tr * zprec_cond(i) / zneb(i))
     404          pfrac_impa(i, kk) = pfrac_impa(i, kk) * (1. - zneb(i) * zfrac_lessi)
     405          frac_impa(i, kk) = 1. - zneb(i) * zfrac_lessi
    409406        END IF
    410407      END DO
    … …  
    420417
    421418  DO i = 1, klon
    422     IF ((t(i,1)+d_t(i,1))<rtt) THEN
     419    IF ((t(i, 1) + d_t(i, 1))<rtt) THEN
    423420      snow(i) = zrfl(i)
    424421    ELSE
    … …  
    427424  END DO
    428425
    429 
    430426END SUBROUTINE fisrtilp_tr
Note: See TracChangeset for help on using the changeset viewer.