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

    r1988 r1992  
    1 !
     1
    22! $Header$
    3 !
    4       SUBROUTINE diagphy(airephy,tit,iprt
    5      $    , tops, topl, sols, soll, sens
    6      $    , evap, rain_fall, snow_fall, ts
    7      $    , d_etp_tot, d_qt_tot, d_ec_tot
    8      $    , fs_bound, fq_bound)
    9 C======================================================================
    10 C
    11 C Purpose:
    12 C    Compute the thermal flux and the watter mass flux at the atmosphere
    13 c    boundaries. Print them and also the atmospheric enthalpy change and
    14 C    the  atmospheric mass change.
    15 C
    16 C Arguments:
    17 C airephy-------input-R-  grid area
    18 C tit---------input-A15- Comment to be added in PRINT (CHARACTER*15)
    19 C iprt--------input-I-  PRINT level ( <=0 : no PRINT)
    20 C tops(klon)--input-R-  SW rad. at TOA (W/m2), positive up.
    21 C topl(klon)--input-R-  LW rad. at TOA (W/m2), positive down
    22 C sols(klon)--input-R-  Net SW flux above surface (W/m2), positive up
    23 C                   (i.e. -1 * flux absorbed by the surface)
    24 C soll(klon)--input-R-  Net LW flux above surface (W/m2), positive up
    25 C                   (i.e. flux emited - flux absorbed by the surface)
    26 C sens(klon)--input-R-  Sensible Flux at surface  (W/m2), positive down
    27 C evap(klon)--input-R-  Evaporation + sublimation watter vapour mass flux
    28 C                   (kg/m2/s), positive up
    29 C rain_fall(klon)
    30 C           --input-R- Liquid  watter mass flux (kg/m2/s), positive down
    31 C snow_fall(klon)
    32 C           --input-R- Solid  watter mass flux (kg/m2/s), positive down
    33 C ts(klon)----input-R- Surface temperature (K)
    34 C d_etp_tot---input-R- Heat flux equivalent to atmospheric enthalpy
    35 C                    change (W/m2)
    36 C d_qt_tot----input-R- Mass flux equivalent to atmospheric watter mass
    37 C                    change (kg/m2/s)
    38 C d_ec_tot----input-R- Flux equivalent to atmospheric cinetic energy
    39 C                    change (W/m2)
    40 C
    41 C fs_bound---output-R- Thermal flux at the atmosphere boundaries (W/m2)
    42 C fq_bound---output-R- Watter mass flux at the atmosphere boundaries (kg/m2/s)
    43 C
    44 C J.L. Dufresne, July 2002
    45 C Version prise sur ~rlmd833/LMDZOR_201102/modipsl/modeles/LMDZ.3.3/libf/phylmd
    46 C  le 25 Novembre 2002.
    47 C======================================================================
    48 C
    49       use dimphy
    50       implicit none
    51 
    52 #include "dimensions.h"
    53 ccccc#include "dimphy.h"
    54 #include "YOMCST.h"
    55 #include "YOETHF.h"
    56 C
    57 C     Input variables
    58       real airephy(klon)
    59       CHARACTER*15 tit
    60       INTEGER iprt
    61       real tops(klon),topl(klon),sols(klon),soll(klon)
    62       real sens(klon),evap(klon),rain_fall(klon),snow_fall(klon)
    63       REAL ts(klon)
    64       REAL d_etp_tot, d_qt_tot, d_ec_tot
    65 c     Output variables
    66       REAL fs_bound, fq_bound
    67 C
    68 C     Local variables
    69       real stops,stopl,ssols,ssoll
    70       real ssens,sfront,slat
    71       real airetot, zcpvap, zcwat, zcice
    72       REAL rain_fall_tot, snow_fall_tot, evap_tot
    73 C
    74       integer i
    75 C
    76       integer pas
    77       save pas
    78       data pas/0/
    79 c$OMP THREADPRIVATE(pas)
    80 C
    81       pas=pas+1
    82       stops=0.
    83       stopl=0.
    84       ssols=0.
    85       ssoll=0.
    86       ssens=0.
    87       sfront = 0.
    88       evap_tot = 0.
    89       rain_fall_tot = 0.
    90       snow_fall_tot = 0.
    91       airetot=0.
    92 C
    93 C     Pour les chaleur specifiques de la vapeur d'eau, de l'eau et de
    94 C     la glace, on travaille par difference a la chaleur specifique de l'
    95 c     air sec. En effet, comme on travaille a niveau de pression donne,
    96 C     toute variation de la masse d'un constituant est totalement
    97 c     compense par une variation de masse d'air.
    98 C
    99       zcpvap=RCPV-RCPD
    100       zcwat=RCW-RCPD
    101       zcice=RCS-RCPD
    102 C
    103       do i=1,klon
    104            stops=stops+tops(i)*airephy(i)
    105            stopl=stopl+topl(i)*airephy(i)
    106            ssols=ssols+sols(i)*airephy(i)
    107            ssoll=ssoll+soll(i)*airephy(i)
    108            ssens=ssens+sens(i)*airephy(i)
    109            sfront = sfront
    110      $         + ( evap(i)*zcpvap-rain_fall(i)*zcwat-snow_fall(i)*zcice
    111      $           ) *ts(i) *airephy(i)
    112            evap_tot = evap_tot + evap(i)*airephy(i)
    113            rain_fall_tot = rain_fall_tot + rain_fall(i)*airephy(i)
    114            snow_fall_tot = snow_fall_tot + snow_fall(i)*airephy(i)
    115            airetot=airetot+airephy(i)
    116       enddo
    117       stops=stops/airetot
    118       stopl=stopl/airetot
    119       ssols=ssols/airetot
    120       ssoll=ssoll/airetot
    121       ssens=ssens/airetot
    122       sfront = sfront/airetot
    123       evap_tot = evap_tot /airetot
    124       rain_fall_tot = rain_fall_tot/airetot
    125       snow_fall_tot = snow_fall_tot/airetot
    126 C
    127       slat = RLVTT * rain_fall_tot + RLSTT * snow_fall_tot
    128 C     Heat flux at atm. boundaries
    129       fs_bound = stops-stopl - (ssols+ssoll)+ssens+sfront
    130      $    + slat
    131 C     Watter flux at atm. boundaries
    132       fq_bound = evap_tot - rain_fall_tot -snow_fall_tot
    133 C
    134       IF (iprt.ge.1) write(6,6666)
    135      $    tit, pas, fs_bound, d_etp_tot, fq_bound, d_qt_tot
    136 C
    137       IF (iprt.ge.1) write(6,6668)
    138      $    tit, pas, d_etp_tot+d_ec_tot-fs_bound, d_qt_tot-fq_bound
    139 C
    140       IF (iprt.ge.2) write(6,6667)
    141      $    tit, pas, stops,stopl,ssols,ssoll,ssens,slat,evap_tot
    142      $    ,rain_fall_tot+snow_fall_tot
    143 
    144       return
    145 
    146  6666 format('Phys. Flux Budget ',a15,1i6,2f8.2,2(1pE13.5))
    147  6667 format('Phys. Boundary Flux ',a15,1i6,6f8.2,2(1pE13.5))
    148  6668 format('Phys. Total Budget ',a15,1i6,f8.2,2(1pE13.5))
    149 
    150       end
    151 
    152 C======================================================================
    153       SUBROUTINE diagetpq(airephy,tit,iprt,idiag,idiag2,dtime
    154      e  ,t,q,ql,qs,u,v,paprs,pplay
    155      s  , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    156 C======================================================================
    157 C
    158 C Purpose:
    159 C    Calcul la difference d'enthalpie et de masse d'eau entre 2 appels,
    160 C    et calcul le flux de chaleur et le flux d'eau necessaire a ces
    161 C    changements. Ces valeurs sont moyennees sur la surface de tout
    162 C    le globe et sont exprime en W/2 et kg/s/m2
    163 C    Outil pour diagnostiquer la conservation de l'energie
    164 C    et de la masse dans la physique. Suppose que les niveau de
    165 c    pression entre couche ne varie pas entre 2 appels.
    166 C
    167 C Plusieurs de ces diagnostics peuvent etre fait en parallele: les
    168 c bilans sont sauvegardes dans des tableaux indices. On parlera
    169 C "d'indice de diagnostic"
    170 c
    171 C
    172 c======================================================================
    173 C Arguments:
    174 C airephy-------input-R-  grid area
    175 C tit-----imput-A15- Comment added in PRINT (CHARACTER*15)
    176 C iprt----input-I-  PRINT level ( <=1 : no PRINT)
    177 C idiag---input-I- indice dans lequel sera range les nouveaux
    178 C                  bilans d' entalpie et de masse
    179 C idiag2--input-I-les nouveaux bilans d'entalpie et de masse
    180 C                 sont compare au bilan de d'enthalpie de masse de
    181 C                 l'indice numero idiag2
    182 C                 Cas parriculier : si idiag2=0, pas de comparaison, on
    183 c                 sort directement les bilans d'enthalpie et de masse
    184 C dtime----input-R- time step (s)
    185 c t--------input-R- temperature (K)
    186 c q--------input-R- vapeur d'eau (kg/kg)
    187 c ql-------input-R- liquid watter (kg/kg)
    188 c qs-------input-R- solid watter (kg/kg)
    189 c u--------input-R- vitesse u
    190 c v--------input-R- vitesse v
    191 c paprs----input-R- pression a intercouche (Pa)
    192 c pplay----input-R- pression au milieu de couche (Pa)
    193 c
    194 C the following total value are computed by UNIT of earth surface
    195 C
    196 C d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
    197 c            change (J/m2) during one time step (dtime) for the whole
    198 C            atmosphere (air, watter vapour, liquid and solid)
    199 C d_qt------output-R- total water mass flux (kg/m2/s) defined as the
    200 C           total watter (kg/m2) change during one time step (dtime),
    201 C d_qw------output-R- same, for the watter vapour only (kg/m2/s)
    202 C d_ql------output-R- same, for the liquid watter only (kg/m2/s)
    203 C d_qs------output-R- same, for the solid watter only (kg/m2/s)
    204 C d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column
    205 C
    206 C     other (COMMON...)
    207 C     RCPD, RCPV, ....
    208 C
    209 C J.L. Dufresne, July 2002
    210 c======================================================================
    211  
    212       USE dimphy
    213       IMPLICIT NONE
    214 C
    215 #include "dimensions.h"
    216 cccccc#include "dimphy.h"
    217 #include "YOMCST.h"
    218 #include "YOETHF.h"
    219 C
    220 c     Input variables
    221       real airephy(klon)
    222       CHARACTER*15 tit
    223       INTEGER iprt,idiag, idiag2
    224       REAL dtime
    225       REAL t(klon,klev), q(klon,klev), ql(klon,klev), qs(klon,klev)
    226       REAL u(klon,klev), v(klon,klev)
    227       REAL paprs(klon,klev+1), pplay(klon,klev)
    228 c     Output variables
    229       REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
    230 C
    231 C     Local variables
    232 c
    233       REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
    234      .  , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
    235 c h_vcol_tot--  total enthalpy of vertical air column
    236 C            (air with watter vapour, liquid and solid) (J/m2)
    237 c h_dair_tot-- total enthalpy of dry air (J/m2)
    238 c h_qw_tot----  total enthalpy of watter vapour (J/m2)
    239 c h_ql_tot----  total enthalpy of liquid watter (J/m2)
    240 c h_qs_tot----  total enthalpy of solid watter  (J/m2)
    241 c qw_tot------  total mass of watter vapour (kg/m2)
    242 c ql_tot------  total mass of liquid watter (kg/m2)
    243 c qs_tot------  total mass of solid watter (kg/m2)
    244 c ec_tot------  total cinetic energy (kg/m2)
    245 C
    246       REAL zairm(klon,klev) ! layer air mass (kg/m2)
    247       REAL  zqw_col(klon)
    248       REAL  zql_col(klon)
    249       REAL  zqs_col(klon)
    250       REAL  zec_col(klon)
    251       REAL  zh_dair_col(klon)
    252       REAL  zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)
    253 C
    254       REAL      d_h_dair, d_h_qw, d_h_ql, d_h_qs
    255 C
    256       REAL airetot, zcpvap, zcwat, zcice
    257 C
    258       INTEGER i, k
    259 C
    260       INTEGER ndiag     ! max number of diagnostic in parallel
    261       PARAMETER (ndiag=10)
    262       integer pas(ndiag)
    263       save pas
    264       data pas/ndiag*0/
    265 c$OMP THREADPRIVATE(pas)
    266 C     
    267       REAL      h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag)
    268      $    , h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag)
    269      $    , ql_pre(ndiag), qs_pre(ndiag) , ec_pre(ndiag)
    270       SAVE      h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre
    271      $        , h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre
    272 c$OMP THREADPRIVATE(h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre)
    273 c$OMP THREADPRIVATE(h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre)
    274 c======================================================================
    275 C
    276       DO k = 1, klev
    277         DO i = 1, klon
    278 C         layer air mass
    279           zairm(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
    280         ENDDO
    281       END DO
    282 C
    283 C     Reset variables
    284       DO i = 1, klon
    285         zqw_col(i)=0.
    286         zql_col(i)=0.
    287         zqs_col(i)=0.
    288         zec_col(i) = 0.
    289         zh_dair_col(i) = 0.
    290         zh_qw_col(i) = 0.
    291         zh_ql_col(i) = 0.
    292         zh_qs_col(i) = 0.
    293       ENDDO
    294 C
    295       zcpvap=RCPV
    296       zcwat=RCW
    297       zcice=RCS
    298 C
    299 C     Compute vertical sum for each atmospheric column
    300 C     ================================================
    301       DO k = 1, klev
    302         DO i = 1, klon
    303 C         Watter mass
    304           zqw_col(i) = zqw_col(i) + q(i,k)*zairm(i,k)
    305           zql_col(i) = zql_col(i) + ql(i,k)*zairm(i,k)
    306           zqs_col(i) = zqs_col(i) + qs(i,k)*zairm(i,k)
    307 C         Cinetic Energy
    308           zec_col(i) =  zec_col(i)
    309      $        +0.5*(u(i,k)**2+v(i,k)**2)*zairm(i,k)
    310 C         Air enthalpy
    311           zh_dair_col(i) = zh_dair_col(i)
    312      $        + RCPD*(1.-q(i,k)-ql(i,k)-qs(i,k))*zairm(i,k)*t(i,k)
    313           zh_qw_col(i) = zh_qw_col(i)
    314      $        + zcpvap*q(i,k)*zairm(i,k)*t(i,k)
    315           zh_ql_col(i) = zh_ql_col(i)
    316      $        + zcwat*ql(i,k)*zairm(i,k)*t(i,k)
    317      $        - RLVTT*ql(i,k)*zairm(i,k)
    318           zh_qs_col(i) = zh_qs_col(i)
    319      $        + zcice*qs(i,k)*zairm(i,k)*t(i,k)
    320      $        - RLSTT*qs(i,k)*zairm(i,k)
    321 
    322         END DO
    323       ENDDO
    324 C
    325 C     Mean over the planete surface
    326 C     =============================
    327       qw_tot = 0.
    328       ql_tot = 0.
    329       qs_tot = 0.
    330       ec_tot = 0.
    331       h_vcol_tot = 0.
    332       h_dair_tot = 0.
    333       h_qw_tot = 0.
    334       h_ql_tot = 0.
    335       h_qs_tot = 0.
    336       airetot=0.
    337 C
    338       do i=1,klon
    339         qw_tot = qw_tot + zqw_col(i)*airephy(i)
    340         ql_tot = ql_tot + zql_col(i)*airephy(i)
    341         qs_tot = qs_tot + zqs_col(i)*airephy(i)
    342         ec_tot = ec_tot + zec_col(i)*airephy(i)
    343         h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
    344         h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
    345         h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
    346         h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
    347         airetot=airetot+airephy(i)
    348       END DO
    349 C
    350       qw_tot = qw_tot/airetot
    351       ql_tot = ql_tot/airetot
    352       qs_tot = qs_tot/airetot
    353       ec_tot = ec_tot/airetot
    354       h_dair_tot = h_dair_tot/airetot
    355       h_qw_tot = h_qw_tot/airetot
    356       h_ql_tot = h_ql_tot/airetot
    357       h_qs_tot = h_qs_tot/airetot
    358 C
    359       h_vcol_tot = h_dair_tot+h_qw_tot+h_ql_tot+h_qs_tot
    360 C
    361 C     Compute the change of the atmospheric state compare to the one
    362 C     stored in "idiag2", and convert it in flux. THis computation
    363 C     is performed IF idiag2 /= 0 and IF it is not the first CALL
    364 c     for "idiag"
    365 C     ===================================
    366 C
    367       IF ( (idiag2.gt.0) .and. (pas(idiag2) .ne. 0) ) THEN
    368         d_h_vcol  = (h_vcol_tot - h_vcol_pre(idiag2) )/dtime
    369         d_h_dair = (h_dair_tot- h_dair_pre(idiag2))/dtime
    370         d_h_qw   = (h_qw_tot  - h_qw_pre(idiag2)  )/dtime
    371         d_h_ql   = (h_ql_tot  - h_ql_pre(idiag2)  )/dtime
    372         d_h_qs   = (h_qs_tot  - h_qs_pre(idiag2)  )/dtime
    373         d_qw     = (qw_tot    - qw_pre(idiag2)    )/dtime
    374         d_ql     = (ql_tot    - ql_pre(idiag2)    )/dtime
    375         d_qs     = (qs_tot    - qs_pre(idiag2)    )/dtime
    376         d_ec     = (ec_tot    - ec_pre(idiag2)    )/dtime
    377         d_qt = d_qw + d_ql + d_qs
    378       ELSE
    379         d_h_vcol = 0.
    380         d_h_dair = 0.
    381         d_h_qw   = 0.
    382         d_h_ql   = 0.
    383         d_h_qs   = 0.
    384         d_qw     = 0.
    385         d_ql     = 0.
    386         d_qs     = 0.
    387         d_ec     = 0.
    388         d_qt     = 0.
    389       ENDIF
    390 C
    391       IF (iprt.ge.2) THEN
    392         WRITE(6,9000) tit,pas(idiag),d_qt,d_qw,d_ql,d_qs
    393  9000   format('Phys. Watter Mass Budget (kg/m2/s)',A15
    394      $      ,1i6,10(1pE14.6))
    395         WRITE(6,9001) tit,pas(idiag), d_h_vcol
    396  9001   format('Phys. Enthalpy Budget (W/m2) ',A15,1i6,10(F8.2))
    397         WRITE(6,9002) tit,pas(idiag), d_ec
    398  9002   format('Phys. Cinetic Energy Budget (W/m2) ',A15,1i6,10(F8.2))
    399       END IF
    400 C
    401 C     Store the new atmospheric state in "idiag"
    402 C
    403       pas(idiag)=pas(idiag)+1
    404       h_vcol_pre(idiag)  = h_vcol_tot
    405       h_dair_pre(idiag) = h_dair_tot
    406       h_qw_pre(idiag)   = h_qw_tot
    407       h_ql_pre(idiag)   = h_ql_tot
    408       h_qs_pre(idiag)   = h_qs_tot
    409       qw_pre(idiag)     = qw_tot
    410       ql_pre(idiag)     = ql_tot
    411       qs_pre(idiag)     = qs_tot
    412       ec_pre (idiag)    = ec_tot
    413 C
    414       RETURN
    415       END
     3
     4SUBROUTINE diagphy(airephy, tit, iprt, tops, topl, sols, soll, sens, evap, &
     5    rain_fall, snow_fall, ts, d_etp_tot, d_qt_tot, d_ec_tot, fs_bound, &
     6    fq_bound)
     7  ! ======================================================================
     8
     9  ! Purpose:
     10  ! Compute the thermal flux and the watter mass flux at the atmosphere
     11  ! boundaries. Print them and also the atmospheric enthalpy change and
     12  ! the  atmospheric mass change.
     13
     14  ! Arguments:
     15  ! airephy-------input-R-  grid area
     16  ! tit---------input-A15- Comment to be added in PRINT (CHARACTER*15)
     17  ! iprt--------input-I-  PRINT level ( <=0 : no PRINT)
     18  ! tops(klon)--input-R-  SW rad. at TOA (W/m2), positive up.
     19  ! topl(klon)--input-R-  LW rad. at TOA (W/m2), positive down
     20  ! sols(klon)--input-R-  Net SW flux above surface (W/m2), positive up
     21  ! (i.e. -1 * flux absorbed by the surface)
     22  ! soll(klon)--input-R-  Net LW flux above surface (W/m2), positive up
     23  ! (i.e. flux emited - flux absorbed by the surface)
     24  ! sens(klon)--input-R-  Sensible Flux at surface  (W/m2), positive down
     25  ! evap(klon)--input-R-  Evaporation + sublimation watter vapour mass flux
     26  ! (kg/m2/s), positive up
     27  ! rain_fall(klon)
     28  ! --input-R- Liquid  watter mass flux (kg/m2/s), positive down
     29  ! snow_fall(klon)
     30  ! --input-R- Solid  watter mass flux (kg/m2/s), positive down
     31  ! ts(klon)----input-R- Surface temperature (K)
     32  ! d_etp_tot---input-R- Heat flux equivalent to atmospheric enthalpy
     33  ! change (W/m2)
     34  ! d_qt_tot----input-R- Mass flux equivalent to atmospheric watter mass
     35  ! change (kg/m2/s)
     36  ! d_ec_tot----input-R- Flux equivalent to atmospheric cinetic energy
     37  ! change (W/m2)
     38
     39  ! fs_bound---output-R- Thermal flux at the atmosphere boundaries (W/m2)
     40  ! fq_bound---output-R- Watter mass flux at the atmosphere boundaries
     41  ! (kg/m2/s)
     42
     43  ! J.L. Dufresne, July 2002
     44  ! Version prise sur
     45  ! ~rlmd833/LMDZOR_201102/modipsl/modeles/LMDZ.3.3/libf/phylmd
     46  ! le 25 Novembre 2002.
     47  ! ======================================================================
     48
     49  USE dimphy
     50  IMPLICIT NONE
     51
     52  include "dimensions.h"
     53  ! cccc#include "dimphy.h"
     54  include "YOMCST.h"
     55  include "YOETHF.h"
     56
     57  ! Input variables
     58  REAL airephy(klon)
     59  CHARACTER *15 tit
     60  INTEGER iprt
     61  REAL tops(klon), topl(klon), sols(klon), soll(klon)
     62  REAL sens(klon), evap(klon), rain_fall(klon), snow_fall(klon)
     63  REAL ts(klon)
     64  REAL d_etp_tot, d_qt_tot, d_ec_tot
     65  ! Output variables
     66  REAL fs_bound, fq_bound
     67
     68  ! Local variables
     69  REAL stops, stopl, ssols, ssoll
     70  REAL ssens, sfront, slat
     71  REAL airetot, zcpvap, zcwat, zcice
     72  REAL rain_fall_tot, snow_fall_tot, evap_tot
     73
     74  INTEGER i
     75
     76  INTEGER pas
     77  SAVE pas
     78  DATA pas/0/
     79  !$OMP THREADPRIVATE(pas)
     80
     81  pas = pas + 1
     82  stops = 0.
     83  stopl = 0.
     84  ssols = 0.
     85  ssoll = 0.
     86  ssens = 0.
     87  sfront = 0.
     88  evap_tot = 0.
     89  rain_fall_tot = 0.
     90  snow_fall_tot = 0.
     91  airetot = 0.
     92
     93  ! Pour les chaleur specifiques de la vapeur d'eau, de l'eau et de
     94  ! la glace, on travaille par difference a la chaleur specifique de l'
     95  ! air sec. En effet, comme on travaille a niveau de pression donne,
     96  ! toute variation de la masse d'un constituant est totalement
     97  ! compense par une variation de masse d'air.
     98
     99  zcpvap = rcpv - rcpd
     100  zcwat = rcw - rcpd
     101  zcice = rcs - rcpd
     102
     103  DO i = 1, klon
     104    stops = stops + tops(i)*airephy(i)
     105    stopl = stopl + topl(i)*airephy(i)
     106    ssols = ssols + sols(i)*airephy(i)
     107    ssoll = ssoll + soll(i)*airephy(i)
     108    ssens = ssens + sens(i)*airephy(i)
     109    sfront = sfront + (evap(i)*zcpvap-rain_fall(i)*zcwat-snow_fall(i)*zcice)* &
     110      ts(i)*airephy(i)
     111    evap_tot = evap_tot + evap(i)*airephy(i)
     112    rain_fall_tot = rain_fall_tot + rain_fall(i)*airephy(i)
     113    snow_fall_tot = snow_fall_tot + snow_fall(i)*airephy(i)
     114    airetot = airetot + airephy(i)
     115  END DO
     116  stops = stops/airetot
     117  stopl = stopl/airetot
     118  ssols = ssols/airetot
     119  ssoll = ssoll/airetot
     120  ssens = ssens/airetot
     121  sfront = sfront/airetot
     122  evap_tot = evap_tot/airetot
     123  rain_fall_tot = rain_fall_tot/airetot
     124  snow_fall_tot = snow_fall_tot/airetot
     125
     126  slat = rlvtt*rain_fall_tot + rlstt*snow_fall_tot
     127  ! Heat flux at atm. boundaries
     128  fs_bound = stops - stopl - (ssols+ssoll) + ssens + sfront + slat
     129  ! Watter flux at atm. boundaries
     130  fq_bound = evap_tot - rain_fall_tot - snow_fall_tot
     131
     132  IF (iprt>=1) WRITE (6, 6666) tit, pas, fs_bound, d_etp_tot, fq_bound, &
     133    d_qt_tot
     134
     135  IF (iprt>=1) WRITE (6, 6668) tit, pas, d_etp_tot + d_ec_tot - fs_bound, &
     136    d_qt_tot - fq_bound
     137
     138  IF (iprt>=2) WRITE (6, 6667) tit, pas, stops, stopl, ssols, ssoll, ssens, &
     139    slat, evap_tot, rain_fall_tot + snow_fall_tot
     140
     141  RETURN
     142
     1436666 FORMAT ('Phys. Flux Budget ', A15, 1I6, 2F8.2, 2(1PE13.5))
     1446667 FORMAT ('Phys. Boundary Flux ', A15, 1I6, 6F8.2, 2(1PE13.5))
     1456668 FORMAT ('Phys. Total Budget ', A15, 1I6, F8.2, 2(1PE13.5))
     146
     147END SUBROUTINE diagphy
     148
     149! ======================================================================
     150SUBROUTINE diagetpq(airephy, tit, iprt, idiag, idiag2, dtime, t, q, ql, qs, &
     151    u, v, paprs, pplay, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     152  ! ======================================================================
     153
     154  ! Purpose:
     155  ! Calcul la difference d'enthalpie et de masse d'eau entre 2 appels,
     156  ! et calcul le flux de chaleur et le flux d'eau necessaire a ces
     157  ! changements. Ces valeurs sont moyennees sur la surface de tout
     158  ! le globe et sont exprime en W/2 et kg/s/m2
     159  ! Outil pour diagnostiquer la conservation de l'energie
     160  ! et de la masse dans la physique. Suppose que les niveau de
     161  ! pression entre couche ne varie pas entre 2 appels.
     162
     163  ! Plusieurs de ces diagnostics peuvent etre fait en parallele: les
     164  ! bilans sont sauvegardes dans des tableaux indices. On parlera
     165  ! "d'indice de diagnostic"
     166
     167
     168  ! ======================================================================
     169  ! Arguments:
     170  ! airephy-------input-R-  grid area
     171  ! tit-----imput-A15- Comment added in PRINT (CHARACTER*15)
     172  ! iprt----input-I-  PRINT level ( <=1 : no PRINT)
     173  ! idiag---input-I- indice dans lequel sera range les nouveaux
     174  ! bilans d' entalpie et de masse
     175  ! idiag2--input-I-les nouveaux bilans d'entalpie et de masse
     176  ! sont compare au bilan de d'enthalpie de masse de
     177  ! l'indice numero idiag2
     178  ! Cas parriculier : si idiag2=0, pas de comparaison, on
     179  ! sort directement les bilans d'enthalpie et de masse
     180  ! dtime----input-R- time step (s)
     181  ! t--------input-R- temperature (K)
     182  ! q--------input-R- vapeur d'eau (kg/kg)
     183  ! ql-------input-R- liquid watter (kg/kg)
     184  ! qs-------input-R- solid watter (kg/kg)
     185  ! u--------input-R- vitesse u
     186  ! v--------input-R- vitesse v
     187  ! paprs----input-R- pression a intercouche (Pa)
     188  ! pplay----input-R- pression au milieu de couche (Pa)
     189
     190  ! the following total value are computed by UNIT of earth surface
     191
     192  ! d_h_vcol--output-R- Heat flux (W/m2) define as the Enthalpy
     193  ! change (J/m2) during one time step (dtime) for the whole
     194  ! atmosphere (air, watter vapour, liquid and solid)
     195  ! d_qt------output-R- total water mass flux (kg/m2/s) defined as the
     196  ! total watter (kg/m2) change during one time step (dtime),
     197  ! d_qw------output-R- same, for the watter vapour only (kg/m2/s)
     198  ! d_ql------output-R- same, for the liquid watter only (kg/m2/s)
     199  ! d_qs------output-R- same, for the solid watter only (kg/m2/s)
     200  ! d_ec------output-R- Cinetic Energy Budget (W/m2) for vertical air column
     201
     202  ! other (COMMON...)
     203  ! RCPD, RCPV, ....
     204
     205  ! J.L. Dufresne, July 2002
     206  ! ======================================================================
     207
     208  USE dimphy
     209  IMPLICIT NONE
     210
     211  include "dimensions.h"
     212  ! ccccc#include "dimphy.h"
     213  include "YOMCST.h"
     214  include "YOETHF.h"
     215
     216  ! Input variables
     217  REAL airephy(klon)
     218  CHARACTER *15 tit
     219  INTEGER iprt, idiag, idiag2
     220  REAL dtime
     221  REAL t(klon, klev), q(klon, klev), ql(klon, klev), qs(klon, klev)
     222  REAL u(klon, klev), v(klon, klev)
     223  REAL paprs(klon, klev+1), pplay(klon, klev)
     224  ! Output variables
     225  REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
     226
     227  ! Local variables
     228
     229  REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot, h_qs_tot, qw_tot, ql_tot, &
     230    qs_tot, ec_tot
     231  ! h_vcol_tot--  total enthalpy of vertical air column
     232  ! (air with watter vapour, liquid and solid) (J/m2)
     233  ! h_dair_tot-- total enthalpy of dry air (J/m2)
     234  ! h_qw_tot----  total enthalpy of watter vapour (J/m2)
     235  ! h_ql_tot----  total enthalpy of liquid watter (J/m2)
     236  ! h_qs_tot----  total enthalpy of solid watter  (J/m2)
     237  ! qw_tot------  total mass of watter vapour (kg/m2)
     238  ! ql_tot------  total mass of liquid watter (kg/m2)
     239  ! qs_tot------  total mass of solid watter (kg/m2)
     240  ! ec_tot------  total cinetic energy (kg/m2)
     241
     242  REAL zairm(klon, klev) ! layer air mass (kg/m2)
     243  REAL zqw_col(klon)
     244  REAL zql_col(klon)
     245  REAL zqs_col(klon)
     246  REAL zec_col(klon)
     247  REAL zh_dair_col(klon)
     248  REAL zh_qw_col(klon), zh_ql_col(klon), zh_qs_col(klon)
     249
     250  REAL d_h_dair, d_h_qw, d_h_ql, d_h_qs
     251
     252  REAL airetot, zcpvap, zcwat, zcice
     253
     254  INTEGER i, k
     255
     256  INTEGER ndiag ! max number of diagnostic in parallel
     257  PARAMETER (ndiag=10)
     258  INTEGER pas(ndiag)
     259  SAVE pas
     260  DATA pas/ndiag*0/
     261  !$OMP THREADPRIVATE(pas)
     262
     263  REAL h_vcol_pre(ndiag), h_dair_pre(ndiag), h_qw_pre(ndiag), &
     264    h_ql_pre(ndiag), h_qs_pre(ndiag), qw_pre(ndiag), ql_pre(ndiag), &
     265    qs_pre(ndiag), ec_pre(ndiag)
     266  SAVE h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre, h_qs_pre, qw_pre, ql_pre, &
     267    qs_pre, ec_pre
     268  !$OMP THREADPRIVATE(h_vcol_pre, h_dair_pre, h_qw_pre, h_ql_pre)
     269  !$OMP THREADPRIVATE(h_qs_pre, qw_pre, ql_pre, qs_pre , ec_pre)
     270  ! ======================================================================
     271
     272  DO k = 1, klev
     273    DO i = 1, klon
     274      ! layer air mass
     275      zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
     276    END DO
     277  END DO
     278
     279  ! Reset variables
     280  DO i = 1, klon
     281    zqw_col(i) = 0.
     282    zql_col(i) = 0.
     283    zqs_col(i) = 0.
     284    zec_col(i) = 0.
     285    zh_dair_col(i) = 0.
     286    zh_qw_col(i) = 0.
     287    zh_ql_col(i) = 0.
     288    zh_qs_col(i) = 0.
     289  END DO
     290
     291  zcpvap = rcpv
     292  zcwat = rcw
     293  zcice = rcs
     294
     295  ! Compute vertical sum for each atmospheric column
     296  ! ================================================
     297  DO k = 1, klev
     298    DO i = 1, klon
     299      ! Watter mass
     300      zqw_col(i) = zqw_col(i) + q(i, k)*zairm(i, k)
     301      zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k)
     302      zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k)
     303      ! Cinetic Energy
     304      zec_col(i) = zec_col(i) + 0.5*(u(i,k)**2+v(i,k)**2)*zairm(i, k)
     305      ! Air enthalpy
     306      zh_dair_col(i) = zh_dair_col(i) + rcpd*(1.-q(i,k)-ql(i,k)-qs(i,k))* &
     307        zairm(i, k)*t(i, k)
     308      zh_qw_col(i) = zh_qw_col(i) + zcpvap*q(i, k)*zairm(i, k)*t(i, k)
     309      zh_ql_col(i) = zh_ql_col(i) + zcwat*ql(i, k)*zairm(i, k)*t(i, k) - &
     310        rlvtt*ql(i, k)*zairm(i, k)
     311      zh_qs_col(i) = zh_qs_col(i) + zcice*qs(i, k)*zairm(i, k)*t(i, k) - &
     312        rlstt*qs(i, k)*zairm(i, k)
     313
     314    END DO
     315  END DO
     316
     317  ! Mean over the planete surface
     318  ! =============================
     319  qw_tot = 0.
     320  ql_tot = 0.
     321  qs_tot = 0.
     322  ec_tot = 0.
     323  h_vcol_tot = 0.
     324  h_dair_tot = 0.
     325  h_qw_tot = 0.
     326  h_ql_tot = 0.
     327  h_qs_tot = 0.
     328  airetot = 0.
     329
     330  DO i = 1, klon
     331    qw_tot = qw_tot + zqw_col(i)*airephy(i)
     332    ql_tot = ql_tot + zql_col(i)*airephy(i)
     333    qs_tot = qs_tot + zqs_col(i)*airephy(i)
     334    ec_tot = ec_tot + zec_col(i)*airephy(i)
     335    h_dair_tot = h_dair_tot + zh_dair_col(i)*airephy(i)
     336    h_qw_tot = h_qw_tot + zh_qw_col(i)*airephy(i)
     337    h_ql_tot = h_ql_tot + zh_ql_col(i)*airephy(i)
     338    h_qs_tot = h_qs_tot + zh_qs_col(i)*airephy(i)
     339    airetot = airetot + airephy(i)
     340  END DO
     341
     342  qw_tot = qw_tot/airetot
     343  ql_tot = ql_tot/airetot
     344  qs_tot = qs_tot/airetot
     345  ec_tot = ec_tot/airetot
     346  h_dair_tot = h_dair_tot/airetot
     347  h_qw_tot = h_qw_tot/airetot
     348  h_ql_tot = h_ql_tot/airetot
     349  h_qs_tot = h_qs_tot/airetot
     350
     351  h_vcol_tot = h_dair_tot + h_qw_tot + h_ql_tot + h_qs_tot
     352
     353  ! Compute the change of the atmospheric state compare to the one
     354  ! stored in "idiag2", and convert it in flux. THis computation
     355  ! is performed IF idiag2 /= 0 and IF it is not the first CALL
     356  ! for "idiag"
     357  ! ===================================
     358
     359  IF ((idiag2>0) .AND. (pas(idiag2)/=0)) THEN
     360    d_h_vcol = (h_vcol_tot-h_vcol_pre(idiag2))/dtime
     361    d_h_dair = (h_dair_tot-h_dair_pre(idiag2))/dtime
     362    d_h_qw = (h_qw_tot-h_qw_pre(idiag2))/dtime
     363    d_h_ql = (h_ql_tot-h_ql_pre(idiag2))/dtime
     364    d_h_qs = (h_qs_tot-h_qs_pre(idiag2))/dtime
     365    d_qw = (qw_tot-qw_pre(idiag2))/dtime
     366    d_ql = (ql_tot-ql_pre(idiag2))/dtime
     367    d_qs = (qs_tot-qs_pre(idiag2))/dtime
     368    d_ec = (ec_tot-ec_pre(idiag2))/dtime
     369    d_qt = d_qw + d_ql + d_qs
     370  ELSE
     371    d_h_vcol = 0.
     372    d_h_dair = 0.
     373    d_h_qw = 0.
     374    d_h_ql = 0.
     375    d_h_qs = 0.
     376    d_qw = 0.
     377    d_ql = 0.
     378    d_qs = 0.
     379    d_ec = 0.
     380    d_qt = 0.
     381  END IF
     382
     383  IF (iprt>=2) THEN
     384    WRITE (6, 9000) tit, pas(idiag), d_qt, d_qw, d_ql, d_qs
     3859000 FORMAT ('Phys. Watter Mass Budget (kg/m2/s)', A15, 1I6, 10(1PE14.6))
     386    WRITE (6, 9001) tit, pas(idiag), d_h_vcol
     3879001 FORMAT ('Phys. Enthalpy Budget (W/m2) ', A15, 1I6, 10(F8.2))
     388    WRITE (6, 9002) tit, pas(idiag), d_ec
     3899002 FORMAT ('Phys. Cinetic Energy Budget (W/m2) ', A15, 1I6, 10(F8.2))
     390  END IF
     391
     392  ! Store the new atmospheric state in "idiag"
     393
     394  pas(idiag) = pas(idiag) + 1
     395  h_vcol_pre(idiag) = h_vcol_tot
     396  h_dair_pre(idiag) = h_dair_tot
     397  h_qw_pre(idiag) = h_qw_tot
     398  h_ql_pre(idiag) = h_ql_tot
     399  h_qs_pre(idiag) = h_qs_tot
     400  qw_pre(idiag) = qw_tot
     401  ql_pre(idiag) = ql_tot
     402  qs_pre(idiag) = qs_tot
     403  ec_pre(idiag) = ec_tot
     404
     405  RETURN
     406END SUBROUTINE diagetpq
Note: See TracChangeset for help on using the changeset viewer.