Ignore:
Timestamp:
Mar 5, 2014, 2:19:12 PM (11 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/tlift.F90

    r1988 r1992  
    1 !
     1
    22! $Header$
    3 !
    4         SUBROUTINE TLIFT(P,T,RR,RS,GZ,PLCL,ICB,NK,
    5      .                  TVP,TPK,CLW,ND,NL,
    6      .                  DTVPDT1,DTVPDQ1)
    7 C
    8 C     Argument NK ajoute (jyg) = Niveau de depart de la
    9 C     convection
    10 C
    11         PARAMETER (NA=60)
    12         REAL GZ(ND),TPK(ND),CLW(ND)
    13         REAL T(ND),RR(ND),RS(ND),TVP(ND),P(ND)
    14         REAL DTVPDT1(ND),DTVPDQ1(ND)   ! Derivatives of parcel virtual
    15 C                                   temperature wrt T1 and Q1
    16 C
    17         REAL CLW_NEW(NA),QI(NA)
    18         REAL DTPDT1(NA),DTPDQ1(NA)      ! Derivatives of parcel temperature
    19 C                                   wrt T1 and Q1
    20  
    21 C
    22         LOGICAL ICE_CONV
    23 C
    24 C   ***   ASSIGN VALUES OF THERMODYNAMIC CONSTANTS     ***
    25 C
    26 c sb        CPD=1005.7
    27 c sb      CPV=1870.0
    28 c sb        CL=4190.0
    29 c sb        CPVMCL=2320.0
    30 c sb        RV=461.5
    31 c sb        RD=287.04
    32 c sb        EPS=RD/RV
    33 c sb        ALV0=2.501E6
    34 ccccccccccccccccccccccc
    35 c constantes coherentes avec le modele du Centre Europeen
    36 c sb      RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644
    37 c sb      RV = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 18.0153
    38 c sb      CPD = 3.5 * RD
    39 c sb      CPV = 4.0 * RV
    40 c sb      CL = 4218.0
    41 c sb      CI=2090.0
    42 c sb      CPVMCL=CL-CPV
    43 c sb      CLMCI=CL-CI
    44 c sb      EPS=RD/RV
    45 c sb      ALV0=2.5008E+06
    46 c sb      ALF0=3.34E+05
    47  
    48 cccccccccccc
    49 c on utilise les constantes thermo du Centre Europeen: (SB)
    50 c
    51 #include "YOMCST.h"
    52        GRAVITY = RG !sb: Pr que gravite ne devienne pas humidite!
    53 c
    54        CPD = RCPD
    55        CPV = RCPV
    56        CL = RCW
    57        CI = RCS
    58        CPVMCL = CL-CPV
    59        CLMCI = CL-CI
    60        EPS = RD/RV
    61        ALV0 = RLVTT
    62        ALF0 = RLMLT ! (ALF0 = RLSTT-RLVTT)
    63 c
    64 cccccccccccccccccccccc
    65 C
    66 C   ***  CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY   ***
    67 C
    68         ICB1=MAX(ICB,2)
    69         ICB1=MIN(ICB,NL)
    70 C
    71 Cjyg1
    72 CC      CPP=CPD*(1.-RR(1))+RR(1)*CPV
    73       CPP=CPD*(1.-RR(NK))+RR(NK)*CPV
    74 Cjyg2
    75       CPINV=1./CPP
    76 Cjyg1
    77 C         ICB may be below condensation level
    78 CCC        DO 100 I=1,ICB1-1
    79 CCC         TPK(I)=T(1)-GZ(I)*CPINV
    80 CCC         TVP(I)=TPK(I)*(1.+RR(1)/EPS)
    81         DO 50 I=1,ICB1
    82          CLW(I)=0.0
    83 50      CONTINUE
    84 C
    85         DO 100 I=NK,ICB1
    86          TPK(I)=T(NK)-(GZ(I) - GZ(NK))*CPINV
    87 Cjyg1
    88 CCC         TVP(I)=TPK(I)*(1.+RR(NK)/EPS)
    89          TVP(I)=TPK(I)*(1.+RR(NK)/EPS-RR(NK))
    90 Cjyg2
    91          DTVPDT1(I) = 1.+RR(NK)/EPS-RR(NK)
    92          DTVPDQ1(I) = TPK(I)*(1./EPS-1.)
    93 C
    94 Cjyg2
    95  
    96   100   CONTINUE
    97  
    98 C
    99 C    ***  FIND LIFTED PARCEL TEMPERATURE AND MIXING RATIO    ***
    100 C
    101 Cjyg1
    102 CC        AH0=(CPD*(1.-RR(1))+CL*RR(1))*T(1)
    103 CC     $     +RR(1)*(ALV0-CPVMCL*(T(1)-273.15))
    104         AH0=(CPD*(1.-RR(NK))+CL*RR(NK))*T(NK)
    105      $     +RR(NK)*(ALV0-CPVMCL*(T(NK)-273.15)) + GZ(NK)
    106 Cjyg2
    107 C
    108 Cjyg1
    109         IMIN = ICB1
    110 C         If ICB is below LCL, start loop at ICB+1
    111         IF (PLCL .LT. P(ICB1)) IMIN = MIN(IMIN+1,NL)
    112 C
    113 CCC        DO 300 I=ICB1,NL
    114         DO 300 I=IMIN,NL
    115 Cjyg2
    116          ALV=ALV0-CPVMCL*(T(I)-273.15)
    117          ALF=ALF0+CLMCI*(T(I)-273.15)
    118  
    119         RG=RS(I)
    120         TG=T(I)
    121 C       S=CPD+ALV*ALV*RG/(RV*T(I)*T(I))
    122 Cjyg1
    123 CC        S=CPD*(1.-RR(1))+CL*RR(1)+ALV*ALV*RG/(RV*T(I)*T(I))
    124         S=CPD*(1.-RR(NK))+CL*RR(NK)+ALV*ALV*RG/(RV*T(I)*T(I))
    125 Cjyg2
    126         S=1./S
    127  
    128         DO 200 J=1,2
    129 Cjyg1
    130 CC         AHG=CPD*TG+(CL-CPD)*RR(1)*TG+ALV*RG+GZ(I)
    131          AHG=CPD*TG+(CL-CPD)*RR(NK)*TG+ALV*RG+GZ(I)
    132 Cjyg2
    133         TG=TG+S*(AH0-AHG)
    134         TC=TG-273.15
    135         DENOM=243.5+TC
    136         DENOM=MAX(DENOM,1.0)
    137 C
    138 C       FORMULE DE BOLTON POUR PSAT
    139 C
    140         ES=6.112*EXP(17.67*TC/DENOM)
    141         RG=EPS*ES/(P(I)-ES*(1.-EPS))
    142  
    143  
    144   200   CONTINUE
    145  
    146 Cjyg1
    147 CC        TPK(I)=(AH0-GZ(I)-ALV*RG)/(CPD+(CL-CPD)*RR(1))
    148         TPK(I)=(AH0-GZ(I)-ALV*RG)/(CPD+(CL-CPD)*RR(NK))
    149 Cjyg2
    150 C       TPK(I)=(AH0-GZ(I)-ALV*RG-(CL-CPD)*T(I)*RR(1))/CPD
    151  
    152 Cjyg1
    153 CC        CLW(I)=RR(1)-RG
    154         CLW(I)=RR(NK)-RG
    155 Cjyg2
    156         CLW(I)=MAX(0.0,CLW(I))
    157 Cjyg1
    158 CCC        TVP(I)=TPK(I)*(1.+RG/EPS)
    159         TVP(I)=TPK(I)*(1.+RG/EPS-RR(NK))
    160 Cjyg2
    161 C
    162 Cjyg1       Derivatives
    163 C
    164         DTPDT1(I) = CPD*S
    165         DTPDQ1(I) = ALV*S
    166 C
    167          DTVPDT1(I) = DTPDT1(I)*(1. + RG/EPS -
    168      .           RR(NK) + ALV*RG/(RD*TPK(I)) )
    169         DTVPDQ1(I) = DTPDQ1(I)*(1. + RG/EPS -
    170      .           RR(NK) + ALV*RG/(RD*TPK(I)) ) - TPK(I)
    171 C
    172 Cjyg2
    173  
    174   300   CONTINUE
    175 C
    176       ICE_CONV = .FALSE.
    177 
    178       IF (ICE_CONV) THEN
    179 C
    180 CJAM
    181 C       RAJOUT DE LA PROCEDURE ICEFRAC
    182 C
    183 c sb        CALL ICEFRAC(T,CLW,CLW_NEW,QI,ND,NL)
    184  
    185         DO 400 I=ICB1,NL
    186         IF (T(I).LT.263.15) THEN
    187         TG=TPK(I)
    188         TC=TPK(I)-273.15
    189         DENOM=243.5+TC
    190         ES=6.112*EXP(17.67*TC/DENOM)
    191         ALV=ALV0-CPVMCL*(T(I)-273.15)
    192         ALF=ALF0+CLMCI*(T(I)-273.15)
    193  
    194         DO J=1,4
    195         ESI=EXP(23.33086-(6111.72784/TPK(I))+0.15215*LOG(TPK(I)))
    196         QSAT_NEW=EPS*ESI/(P(I)-ESI*(1.-EPS))
    197 CCC        SNEW= CPD*(1.-RR(1))+CL*RR(1)+ALV*ALV*QSAT_NEW/(RV*TPK(I)*TPK(I))
    198         SNEW= CPD*(1.-RR(NK))+CL*RR(NK)
    199      .        +ALV*ALV*QSAT_NEW/(RV*TPK(I)*TPK(I))
    200 C
    201         SNEW=1./SNEW
    202         TPK(I)=TG+(ALF*QI(I)+ALV*RG*(1.-(ESI/ES)))*SNEW
    203 c@$$        PRINT*,'################################'
    204 c@$$        PRINT*,TPK(I)
    205 c@$$        PRINT*,(ALF*QI(I)+ALV*RG*(1.-(ESI/ES)))*SNEW
    206         ENDDO
    207 CCC        CLW(I)=RR(1)-QSAT_NEW
    208         CLW(I)=RR(NK)-QSAT_NEW
    209         CLW(I)=MAX(0.0,CLW(I))
    210 Cjyg1
    211 CCC        TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS)
    212         TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS-RR(NK))
    213 Cjyg2
    214         ELSE
     3
     4SUBROUTINE tlift(p, t, rr, rs, gz, plcl, icb, nk, tvp, tpk, clw, nd, nl, &
     5    dtvpdt1, dtvpdq1)
     6
     7  ! Argument NK ajoute (jyg) = Niveau de depart de la
     8  ! convection
     9
     10  PARAMETER (na=60)
     11  REAL gz(nd), tpk(nd), clw(nd)
     12  REAL t(nd), rr(nd), rs(nd), tvp(nd), p(nd)
     13  REAL dtvpdt1(nd), dtvpdq1(nd) ! Derivatives of parcel virtual
     14  ! temperature wrt T1 and Q1
     15
     16  REAL clw_new(na), qi(na)
     17  REAL dtpdt1(na), dtpdq1(na) ! Derivatives of parcel temperature
     18  ! wrt T1 and Q1
     19
     20
     21  LOGICAL ice_conv
     22
     23  ! ***   ASSIGN VALUES OF THERMODYNAMIC CONSTANTS     ***
     24
     25  ! sb        CPD=1005.7
     26  ! sb      CPV=1870.0
     27  ! sb        CL=4190.0
     28  ! sb        CPVMCL=2320.0
     29  ! sb        RV=461.5
     30  ! sb        RD=287.04
     31  ! sb        EPS=RD/RV
     32  ! sb        ALV0=2.501E6
     33  ! cccccccccccccccccccccc
     34  ! constantes coherentes avec le modele du Centre Europeen
     35  ! sb      RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644
     36  ! sb      RV = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 18.0153
     37  ! sb      CPD = 3.5 * RD
     38  ! sb      CPV = 4.0 * RV
     39  ! sb      CL = 4218.0
     40  ! sb      CI=2090.0
     41  ! sb      CPVMCL=CL-CPV
     42  ! sb      CLMCI=CL-CI
     43  ! sb      EPS=RD/RV
     44  ! sb      ALV0=2.5008E+06
     45  ! sb      ALF0=3.34E+05
     46
     47  ! ccccccccccc
     48  ! on utilise les constantes thermo du Centre Europeen: (SB)
     49
     50  include "YOMCST.h"
     51  gravity = rg !sb: Pr que gravite ne devienne pas humidite!
     52
     53  cpd = rcpd
     54  cpv = rcpv
     55  cl = rcw
     56  ci = rcs
     57  cpvmcl = cl - cpv
     58  clmci = cl - ci
     59  eps = rd/rv
     60  alv0 = rlvtt
     61  alf0 = rlmlt ! (ALF0 = RLSTT-RLVTT)
     62
     63  ! ccccccccccccccccccccc
     64
     65  ! ***  CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY   ***
     66
     67  icb1 = max(icb, 2)
     68  icb1 = min(icb, nl)
     69
     70  ! jyg1
     71  ! C      CPP=CPD*(1.-RR(1))+RR(1)*CPV
     72  cpp = cpd*(1.-rr(nk)) + rr(nk)*cpv
     73  ! jyg2
     74  cpinv = 1./cpp
     75  ! jyg1
     76  ! ICB may be below condensation level
     77  ! CC        DO 100 I=1,ICB1-1
     78  ! CC         TPK(I)=T(1)-GZ(I)*CPINV
     79  ! CC         TVP(I)=TPK(I)*(1.+RR(1)/EPS)
     80  DO i = 1, icb1
     81    clw(i) = 0.0
     82  END DO
     83
     84  DO i = nk, icb1
     85    tpk(i) = t(nk) - (gz(i)-gz(nk))*cpinv
     86    ! jyg1
     87    ! CC         TVP(I)=TPK(I)*(1.+RR(NK)/EPS)
     88    tvp(i) = tpk(i)*(1.+rr(nk)/eps-rr(nk))
     89    ! jyg2
     90    dtvpdt1(i) = 1. + rr(nk)/eps - rr(nk)
     91    dtvpdq1(i) = tpk(i)*(1./eps-1.)
     92
     93    ! jyg2
     94
     95  END DO
     96
     97
     98  ! ***  FIND LIFTED PARCEL TEMPERATURE AND MIXING RATIO    ***
     99
     100  ! jyg1
     101  ! C        AH0=(CPD*(1.-RR(1))+CL*RR(1))*T(1)
     102  ! C     $     +RR(1)*(ALV0-CPVMCL*(T(1)-273.15))
     103  ah0 = (cpd*(1.-rr(nk))+cl*rr(nk))*t(nk) + rr(nk)*(alv0-cpvmcl*(t(nk)-273.15 &
     104    )) + gz(nk)
     105  ! jyg2
     106
     107  ! jyg1
     108  imin = icb1
     109  ! If ICB is below LCL, start loop at ICB+1
     110  IF (plcl<p(icb1)) imin = min(imin+1, nl)
     111
     112  ! CC        DO 300 I=ICB1,NL
     113  DO i = imin, nl
     114    ! jyg2
     115    alv = alv0 - cpvmcl*(t(i)-273.15)
     116    alf = alf0 + clmci*(t(i)-273.15)
     117
     118    rg = rs(i)
     119    tg = t(i)
     120    ! S=CPD+ALV*ALV*RG/(RV*T(I)*T(I))
     121    ! jyg1
     122    ! C        S=CPD*(1.-RR(1))+CL*RR(1)+ALV*ALV*RG/(RV*T(I)*T(I))
     123    s = cpd*(1.-rr(nk)) + cl*rr(nk) + alv*alv*rg/(rv*t(i)*t(i))
     124    ! jyg2
     125    s = 1./s
     126
     127    DO j = 1, 2
     128      ! jyg1
     129      ! C         AHG=CPD*TG+(CL-CPD)*RR(1)*TG+ALV*RG+GZ(I)
     130      ahg = cpd*tg + (cl-cpd)*rr(nk)*tg + alv*rg + gz(i)
     131      ! jyg2
     132      tg = tg + s*(ah0-ahg)
     133      tc = tg - 273.15
     134      denom = 243.5 + tc
     135      denom = max(denom, 1.0)
     136
     137      ! FORMULE DE BOLTON POUR PSAT
     138
     139      es = 6.112*exp(17.67*tc/denom)
     140      rg = eps*es/(p(i)-es*(1.-eps))
     141
     142
     143    END DO
     144
     145    ! jyg1
     146    ! C        TPK(I)=(AH0-GZ(I)-ALV*RG)/(CPD+(CL-CPD)*RR(1))
     147    tpk(i) = (ah0-gz(i)-alv*rg)/(cpd+(cl-cpd)*rr(nk))
     148    ! jyg2
     149    ! TPK(I)=(AH0-GZ(I)-ALV*RG-(CL-CPD)*T(I)*RR(1))/CPD
     150
     151    ! jyg1
     152    ! C        CLW(I)=RR(1)-RG
     153    clw(i) = rr(nk) - rg
     154    ! jyg2
     155    clw(i) = max(0.0, clw(i))
     156    ! jyg1
     157    ! CC        TVP(I)=TPK(I)*(1.+RG/EPS)
     158    tvp(i) = tpk(i)*(1.+rg/eps-rr(nk))
     159    ! jyg2
     160
     161    ! jyg1       Derivatives
     162
     163    dtpdt1(i) = cpd*s
     164    dtpdq1(i) = alv*s
     165
     166    dtvpdt1(i) = dtpdt1(i)*(1.+rg/eps-rr(nk)+alv*rg/(rd*tpk(i)))
     167    dtvpdq1(i) = dtpdq1(i)*(1.+rg/eps-rr(nk)+alv*rg/(rd*tpk(i))) - tpk(i)
     168
     169    ! jyg2
     170
     171  END DO
     172
     173  ice_conv = .FALSE.
     174
     175  IF (ice_conv) THEN
     176
     177    ! JAM
     178    ! RAJOUT DE LA PROCEDURE ICEFRAC
     179
     180    ! sb        CALL ICEFRAC(T,CLW,CLW_NEW,QI,ND,NL)
     181
     182    DO i = icb1, nl
     183      IF (t(i)<263.15) THEN
     184        tg = tpk(i)
     185        tc = tpk(i) - 273.15
     186        denom = 243.5 + tc
     187        es = 6.112*exp(17.67*tc/denom)
     188        alv = alv0 - cpvmcl*(t(i)-273.15)
     189        alf = alf0 + clmci*(t(i)-273.15)
     190
     191        DO j = 1, 4
     192          esi = exp(23.33086-(6111.72784/tpk(i))+0.15215*log(tpk(i)))
     193          qsat_new = eps*esi/(p(i)-esi*(1.-eps))
     194          ! CC        SNEW=
     195          ! CPD*(1.-RR(1))+CL*RR(1)+ALV*ALV*QSAT_NEW/(RV*TPK(I)*TPK(I))
     196          snew = cpd*(1.-rr(nk)) + cl*rr(nk) + alv*alv*qsat_new/(rv*tpk(i)* &
     197            tpk(i))
     198
     199          snew = 1./snew
     200          tpk(i) = tg + (alf*qi(i)+alv*rg*(1.-(esi/es)))*snew
     201          ! @$$        PRINT*,'################################'
     202          ! @$$        PRINT*,TPK(I)
     203          ! @$$        PRINT*,(ALF*QI(I)+ALV*RG*(1.-(ESI/ES)))*SNEW
     204        END DO
     205        ! CC        CLW(I)=RR(1)-QSAT_NEW
     206        clw(i) = rr(nk) - qsat_new
     207        clw(i) = max(0.0, clw(i))
     208        ! jyg1
     209        ! CC        TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS)
     210        tvp(i) = tpk(i)*(1.+qsat_new/eps-rr(nk))
     211        ! jyg2
     212      ELSE
    215213        CONTINUE
    216         ENDIF
    217  
    218   400   CONTINUE
    219 C
    220       ENDIF
    221 C
    222  
    223 ******************************************************
    224 ** BK :  RAJOUT DE LA TEMPERATURE DES ASCENDANCES
    225 **   NON DILUES AU  NIVEAU KLEV = ND
    226 **   POSONS LE ENVIRON EGAL A CELUI DE KLEV-1
    227 ********************************************************
    228  
    229       TPK(NL+1)=TPK(NL)
    230  
    231 *******************************************************
    232 
    233       RG = GRAVITY ! RG redevient la gravite de YOMCST (sb)
    234  
    235  
    236         RETURN
    237         END
    238  
    239  
    240  
    241  
    242  
    243  
    244  
    245  
     214      END IF
     215
     216    END DO
     217
     218  END IF
     219
     220
     221  ! *****************************************************
     222  ! * BK :  RAJOUT DE LA TEMPERATURE DES ASCENDANCES
     223  ! *   NON DILUES AU  NIVEAU KLEV = ND
     224  ! *   POSONS LE ENVIRON EGAL A CELUI DE KLEV-1
     225  ! *******************************************************
     226
     227  tpk(nl+1) = tpk(nl)
     228
     229  ! ******************************************************
     230
     231  rg = gravity ! RG redevient la gravite de YOMCST (sb)
     232
     233
     234  RETURN
     235END SUBROUTINE tlift
     236
     237
     238
     239
     240
     241
     242
     243
Note: See TracChangeset for help on using the changeset viewer.