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

    r1988 r1992  
    1 !
     1
    22! $Id$
    3 !
    4       SUBROUTINE conema3 (dtime,paprs,pplay,t,q,u,v,tra,ntra,
    5      .             work1,work2,d_t,d_q,d_u,d_v,d_tra,
    6      .             rain, snow, kbas, ktop,
    7      .             upwd,dnwd,dnwdbis,bas,top,Ma,cape,tvp,rflag,
    8      .             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
    9      .             qcond_incld)
    10 
    11       USE dimphy
    12       USE infotrac, ONLY : nbtr
    13       IMPLICIT none
    14 c======================================================================
    15 c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
    16 c Objet: schema de convection de Emanuel (1991) interface
    17 c Mai 1998: Interface modifiee pour implementation dans LMDZ
    18 c======================================================================
    19 c Arguments:
    20 c dtime---input-R-pas d'integration (s)
    21 c paprs---input-R-pression inter-couches (Pa)
    22 c pplay---input-R-pression au milieu des couches (Pa)
    23 c t-------input-R-temperature (K)
    24 c q-------input-R-humidite specifique (kg/kg)
    25 c u-------input-R-vitesse du vent zonal (m/s)
    26 c v-------input-R-vitesse duvent meridien (m/s)
    27 c tra-----input-R-tableau de rapport de melange des traceurs
    28 c work*: input et output: deux variables de travail,
    29 c                            on peut les mettre a 0 au debut
    30 c
    31 C d_t-----output-R-increment de la temperature
    32 c d_q-----output-R-increment de la vapeur d'eau
    33 c d_u-----output-R-increment de la vitesse zonale
    34 c d_v-----output-R-increment de la vitesse meridienne
    35 c d_tra---output-R-increment du contenu en traceurs
    36 c rain----output-R-la pluie (mm/s)
    37 c snow----output-R-la neige (mm/s)
    38 c kbas----output-R-bas du nuage (integer)
    39 c ktop----output-R-haut du nuage (integer)
    40 c upwd----output-R-saturated updraft mass flux (kg/m**2/s)
    41 c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
    42 c dnwdbis-output-R-unsaturated downdraft mass flux (kg/m**2/s)
    43 c bas-----output-R-bas du nuage (real)
    44 c top-----output-R-haut du nuage (real)
    45 c Ma------output-R-flux ascendant non dilue (kg/m**2/s)
    46 c cape----output-R-CAPE
    47 c tvp-----output-R-virtual temperature of the lifted parcel
    48 c rflag---output-R-flag sur le fonctionnement de convect
    49 c pbase---output-R-pression a la base du nuage (Pa)
    50 c bbase---output-R-buoyancy a la base du nuage (K)
    51 c dtvpdt1-output-R-derivative of parcel virtual temp wrt T1
    52 c dtvpdq1-output-R-derivative of parcel virtual temp wrt Q1
    53 c dplcldt-output-R-derivative of the PCP pressure wrt T1
    54 c dplcldr-output-R-derivative of the PCP pressure wrt Q1
    55 c======================================================================
    56 c
    57 #include "dimensions.h"
    58 #include "conema3.h"
    59       INTEGER i, l,m,itra
    60       INTEGER ntra       ! if no tracer transport
    61                          ! is needed, set ntra = 1 (or 0)
    62       REAL dtime
    63 c
    64       REAL d_t2(klon,klev), d_q2(klon,klev) ! sbl
    65       REAL d_u2(klon,klev), d_v2(klon,klev) ! sbl
    66       REAL em_d_t2(klev), em_d_q2(klev)     ! sbl   
    67       REAL em_d_u2(klev), em_d_v2(klev)     ! sbl   
    68 c
    69       REAL paprs(klon,klev+1), pplay(klon,klev)
    70       REAL t(klon,klev), q(klon,klev), d_t(klon,klev), d_q(klon,klev)
    71       REAL u(klon,klev), v(klon,klev), tra(klon,klev,ntra)
    72       REAL d_u(klon,klev), d_v(klon,klev), d_tra(klon,klev,ntra)
    73       REAL work1(klon,klev), work2(klon,klev)
    74       REAL upwd(klon,klev), dnwd(klon,klev), dnwdbis(klon,klev)
    75       REAL rain(klon)
    76       REAL snow(klon)
    77       REAL cape(klon), tvp(klon,klev), rflag(klon)
    78       REAL pbase(klon), bbase(klon)
    79       REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
    80       REAL dplcldt(klon), dplcldr(klon)
    81       INTEGER kbas(klon), ktop(klon)
    82 
    83       REAL wd(klon)
    84       REAL qcond_incld(klon,klev)
    85 c
    86       LOGICAL,SAVE :: first=.true.
    87 c$OMP THREADPRIVATE(first)
    88      
    89 cym      REAL em_t(klev)
    90       REAL,ALLOCATABLE,SAVE :: em_t(:)
    91 c$OMP THREADPRIVATE(em_t) 
    92 cym      REAL em_q(klev)
    93       REAL,ALLOCATABLE,SAVE :: em_q(:)
    94 c$OMP THREADPRIVATE(em_q)
    95 cym      REAL em_qs(klev)
    96       REAL,ALLOCATABLE,SAVE :: em_qs(:)
    97 c$OMP THREADPRIVATE(em_qs) 
    98 cym      REAL em_u(klev), em_v(klev), em_tra(klev,nbtr)
    99       REAL,ALLOCATABLE,SAVE :: em_u(:),em_v(:),em_tra(:,:)
    100 c$OMP THREADPRIVATE(em_u,em_v,em_tra)     
    101 cym      REAL em_ph(klev+1), em_p(klev)
    102       REAL,ALLOCATABLE,SAVE ::em_ph(:),em_p(:)
    103 c$OMP THREADPRIVATE(em_ph,em_p)
    104 cym      REAL em_work1(klev), em_work2(klev)
    105       REAL,ALLOCATABLE,SAVE ::em_work1(:),em_work2(:)
    106 c$OMP THREADPRIVATE(em_work1,em_work2)     
    107 cym      REAL em_precip, em_d_t(klev), em_d_q(klev)
    108       REAL,SAVE :: em_precip
    109 c$OMP THREADPRIVATE(em_precip)     
    110       REAL,ALLOCATABLE,SAVE :: em_d_t(:),em_d_q(:)
    111 c$OMP THREADPRIVATE(em_d_t,em_d_q)
    112 cym      REAL em_d_u(klev), em_d_v(klev), em_d_tra(klev,nbtr)
    113       REAL,ALLOCATABLE,SAVE ::em_d_u(:),em_d_v(:),em_d_tra(:,:)
    114 c$OMP THREADPRIVATE(em_d_u,em_d_v,em_d_tra)     
    115 cym      REAL em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev)
    116       REAL,ALLOCATABLE,SAVE :: em_upwd(:),em_dnwd(:),em_dnwdbis(:)
    117 c$OMP THREADPRIVATE(em_upwd,em_dnwd,em_dnwdbis)
    118       REAL em_dtvpdt1(klev), em_dtvpdq1(klev)
    119       REAL em_dplcldt, em_dplcldr
    120 cym      SAVE em_t,em_q, em_qs, em_ph, em_p, em_work1, em_work2
    121 cym      SAVE em_u,em_v, em_tra
    122 cym      SAVE em_d_u,em_d_v, em_d_tra
    123 cym      SAVE em_precip, em_d_t, em_d_q, em_upwd, em_dnwd, em_dnwdbis
    124 
    125       INTEGER em_bas, em_top
    126       SAVE em_bas, em_top
    127 c$OMP THREADPRIVATE(em_bas,em_top)
    128       REAL em_wd
    129       REAL em_qcond(klev)
    130       REAL em_qcondc(klev)
    131 c
    132       REAL zx_t, zx_qs, zdelta, zcor
    133       INTEGER iflag
    134       REAL sigsum
    135 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    136 c     VARIABLES A SORTIR
    137 cccccccccccccccccccccccccccccccccccccccccccccccccc
    138  
    139 cym      REAL emmip(klev) !variation de flux ascnon dilue i et i+1
    140       REAL,ALLOCATABLE,SAVE ::emmip(:)
    141 c$OMP THREADPRIVATE(emmip)
    142 cym      SAVE emmip
    143 cym      real emMke(klev)
    144       REAL,ALLOCATABLE,SAVE ::emMke(:)
    145 c$OMP THREADPRIVATE(emMke)
    146 cym      save emMke
    147       real top
    148       real bas
    149 cym      real emMa(klev)
    150       REAL,ALLOCATABLE,SAVE ::emMa(:)
    151 c$OMP THREADPRIVATE(emMa)
    152 cym      save emMa
    153       real Ma(klon,klev)
    154       real Ment(klev,klev)
    155       real Qent(klev,klev)
    156       real TPS(klev),TLS(klev)
    157       real SIJ(klev,klev)
    158       real em_CAPE, em_TVP(klev)
    159       real em_pbase, em_bbase
    160       integer iw,j,k,ix,iy
    161 
    162 c -- sb: pour schema nuages:
    163 
    164        integer iflagcon
    165        integer em_ifc(klev)
    166      
    167        real em_pradj
    168        real em_cldf(klev), em_cldq(klev)
    169        real em_ftadj(klev), em_fradj(klev)
    170 
    171        integer ifc(klon,klev)
    172        real pradj(klon)
    173        real cldf(klon,klev), cldq(klon,klev)
    174        real ftadj(klon,klev), fqadj(klon,klev)
    175 
    176 c sb --
    177 
    178 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    179 c
    180 #include "YOMCST.h"
    181 #include "YOETHF.h"
    182 #include "FCTTRE.h"
    183      
    184       if (first) then
    185  
    186         allocate(em_t(klev))
    187         allocate(em_q(klev))
    188         allocate(em_qs(klev))
    189         allocate(em_u(klev), em_v(klev), em_tra(klev,nbtr))
    190         allocate(em_ph(klev+1), em_p(klev))
    191         allocate(em_work1(klev), em_work2(klev))
    192         allocate(em_d_t(klev), em_d_q(klev))
    193         allocate(em_d_u(klev), em_d_v(klev), em_d_tra(klev,nbtr))
    194         allocate(em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev))
    195         allocate(emmip(klev))
    196         allocate(emMke(klev))
    197         allocate(emMa(klev))
    198  
    199         first=.false.
    200       endif
    201  
    202       qcond_incld(:,:) = 0.
    203 c
    204 c@$$      print*,'debut conema'
    205 
    206       DO 999 i = 1, klon
    207       DO l = 1, klev+1
    208          em_ph(l) = paprs(i,l) / 100.0
    209       ENDDO
    210 c
     3
     4SUBROUTINE conema3(dtime, paprs, pplay, t, q, u, v, tra, ntra, work1, work2, &
     5    d_t, d_q, d_u, d_v, d_tra, rain, snow, kbas, ktop, upwd, dnwd, dnwdbis, &
     6    bas, top, ma, cape, tvp, rflag, pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, &
     7    dplcldr, qcond_incld)
     8
     9  USE dimphy
     10  USE infotrac, ONLY: nbtr
     11  IMPLICIT NONE
     12  ! ======================================================================
     13  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
     14  ! Objet: schema de convection de Emanuel (1991) interface
     15  ! Mai 1998: Interface modifiee pour implementation dans LMDZ
     16  ! ======================================================================
     17  ! Arguments:
     18  ! dtime---input-R-pas d'integration (s)
     19  ! paprs---input-R-pression inter-couches (Pa)
     20  ! pplay---input-R-pression au milieu des couches (Pa)
     21  ! t-------input-R-temperature (K)
     22  ! q-------input-R-humidite specifique (kg/kg)
     23  ! u-------input-R-vitesse du vent zonal (m/s)
     24  ! v-------input-R-vitesse duvent meridien (m/s)
     25  ! tra-----input-R-tableau de rapport de melange des traceurs
     26  ! work*: input et output: deux variables de travail,
     27  ! on peut les mettre a 0 au debut
     28
     29  ! d_t-----output-R-increment de la temperature
     30  ! d_q-----output-R-increment de la vapeur d'eau
     31  ! d_u-----output-R-increment de la vitesse zonale
     32  ! d_v-----output-R-increment de la vitesse meridienne
     33  ! d_tra---output-R-increment du contenu en traceurs
     34  ! rain----output-R-la pluie (mm/s)
     35  ! snow----output-R-la neige (mm/s)
     36  ! kbas----output-R-bas du nuage (integer)
     37  ! ktop----output-R-haut du nuage (integer)
     38  ! upwd----output-R-saturated updraft mass flux (kg/m**2/s)
     39  ! dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
     40  ! dnwdbis-output-R-unsaturated downdraft mass flux (kg/m**2/s)
     41  ! bas-----output-R-bas du nuage (real)
     42  ! top-----output-R-haut du nuage (real)
     43  ! Ma------output-R-flux ascendant non dilue (kg/m**2/s)
     44  ! cape----output-R-CAPE
     45  ! tvp-----output-R-virtual temperature of the lifted parcel
     46  ! rflag---output-R-flag sur le fonctionnement de convect
     47  ! pbase---output-R-pression a la base du nuage (Pa)
     48  ! bbase---output-R-buoyancy a la base du nuage (K)
     49  ! dtvpdt1-output-R-derivative of parcel virtual temp wrt T1
     50  ! dtvpdq1-output-R-derivative of parcel virtual temp wrt Q1
     51  ! dplcldt-output-R-derivative of the PCP pressure wrt T1
     52  ! dplcldr-output-R-derivative of the PCP pressure wrt Q1
     53  ! ======================================================================
     54
     55  include "dimensions.h"
     56  include "conema3.h"
     57  INTEGER i, l, m, itra
     58  INTEGER ntra ! if no tracer transport
     59    ! is needed, set ntra = 1 (or 0)
     60  REAL dtime
     61
     62  REAL d_t2(klon, klev), d_q2(klon, klev) ! sbl
     63  REAL d_u2(klon, klev), d_v2(klon, klev) ! sbl
     64  REAL em_d_t2(klev), em_d_q2(klev) ! sbl
     65  REAL em_d_u2(klev), em_d_v2(klev) ! sbl
     66
     67  REAL paprs(klon, klev+1), pplay(klon, klev)
     68  REAL t(klon, klev), q(klon, klev), d_t(klon, klev), d_q(klon, klev)
     69  REAL u(klon, klev), v(klon, klev), tra(klon, klev, ntra)
     70  REAL d_u(klon, klev), d_v(klon, klev), d_tra(klon, klev, ntra)
     71  REAL work1(klon, klev), work2(klon, klev)
     72  REAL upwd(klon, klev), dnwd(klon, klev), dnwdbis(klon, klev)
     73  REAL rain(klon)
     74  REAL snow(klon)
     75  REAL cape(klon), tvp(klon, klev), rflag(klon)
     76  REAL pbase(klon), bbase(klon)
     77  REAL dtvpdt1(klon, klev), dtvpdq1(klon, klev)
     78  REAL dplcldt(klon), dplcldr(klon)
     79  INTEGER kbas(klon), ktop(klon)
     80
     81  REAL wd(klon)
     82  REAL qcond_incld(klon, klev)
     83
     84  LOGICAL, SAVE :: first = .TRUE.
     85  !$OMP THREADPRIVATE(first)
     86
     87  ! ym      REAL em_t(klev)
     88  REAL, ALLOCATABLE, SAVE :: em_t(:)
     89  !$OMP THREADPRIVATE(em_t)
     90  ! ym      REAL em_q(klev)
     91  REAL, ALLOCATABLE, SAVE :: em_q(:)
     92  !$OMP THREADPRIVATE(em_q)
     93  ! ym      REAL em_qs(klev)
     94  REAL, ALLOCATABLE, SAVE :: em_qs(:)
     95  !$OMP THREADPRIVATE(em_qs)
     96  ! ym      REAL em_u(klev), em_v(klev), em_tra(klev,nbtr)
     97  REAL, ALLOCATABLE, SAVE :: em_u(:), em_v(:), em_tra(:, :)
     98  !$OMP THREADPRIVATE(em_u,em_v,em_tra)
     99  ! ym      REAL em_ph(klev+1), em_p(klev)
     100  REAL, ALLOCATABLE, SAVE :: em_ph(:), em_p(:)
     101  !$OMP THREADPRIVATE(em_ph,em_p)
     102  ! ym      REAL em_work1(klev), em_work2(klev)
     103  REAL, ALLOCATABLE, SAVE :: em_work1(:), em_work2(:)
     104  !$OMP THREADPRIVATE(em_work1,em_work2)
     105  ! ym      REAL em_precip, em_d_t(klev), em_d_q(klev)
     106  REAL, SAVE :: em_precip
     107  !$OMP THREADPRIVATE(em_precip)
     108  REAL, ALLOCATABLE, SAVE :: em_d_t(:), em_d_q(:)
     109  !$OMP THREADPRIVATE(em_d_t,em_d_q)
     110  ! ym      REAL em_d_u(klev), em_d_v(klev), em_d_tra(klev,nbtr)
     111  REAL, ALLOCATABLE, SAVE :: em_d_u(:), em_d_v(:), em_d_tra(:, :)
     112  !$OMP THREADPRIVATE(em_d_u,em_d_v,em_d_tra)
     113  ! ym      REAL em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev)
     114  REAL, ALLOCATABLE, SAVE :: em_upwd(:), em_dnwd(:), em_dnwdbis(:)
     115  !$OMP THREADPRIVATE(em_upwd,em_dnwd,em_dnwdbis)
     116  REAL em_dtvpdt1(klev), em_dtvpdq1(klev)
     117  REAL em_dplcldt, em_dplcldr
     118  ! ym      SAVE em_t,em_q, em_qs, em_ph, em_p, em_work1, em_work2
     119  ! ym      SAVE em_u,em_v, em_tra
     120  ! ym      SAVE em_d_u,em_d_v, em_d_tra
     121  ! ym      SAVE em_precip, em_d_t, em_d_q, em_upwd, em_dnwd, em_dnwdbis
     122
     123  INTEGER em_bas, em_top
     124  SAVE em_bas, em_top
     125  !$OMP THREADPRIVATE(em_bas,em_top)
     126  REAL em_wd
     127  REAL em_qcond(klev)
     128  REAL em_qcondc(klev)
     129
     130  REAL zx_t, zx_qs, zdelta, zcor
     131  INTEGER iflag
     132  REAL sigsum
     133  ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     134  ! VARIABLES A SORTIR
     135  ! ccccccccccccccccccccccccccccccccccccccccccccccccc
     136
     137  ! ym      REAL emmip(klev) !variation de flux ascnon dilue i et i+1
     138  REAL, ALLOCATABLE, SAVE :: emmip(:)
     139  !$OMP THREADPRIVATE(emmip)
     140  ! ym      SAVE emmip
     141  ! ym      real emMke(klev)
     142  REAL, ALLOCATABLE, SAVE :: emmke(:)
     143  !$OMP THREADPRIVATE(emMke)
     144  ! ym      save emMke
     145  REAL top
     146  REAL bas
     147  ! ym      real emMa(klev)
     148  REAL, ALLOCATABLE, SAVE :: emma(:)
     149  !$OMP THREADPRIVATE(emMa)
     150  ! ym      save emMa
     151  REAL ma(klon, klev)
     152  REAL ment(klev, klev)
     153  REAL qent(klev, klev)
     154  REAL tps(klev), tls(klev)
     155  REAL sij(klev, klev)
     156  REAL em_cape, em_tvp(klev)
     157  REAL em_pbase, em_bbase
     158  INTEGER iw, j, k, ix, iy
     159
     160  ! -- sb: pour schema nuages:
     161
     162  INTEGER iflagcon
     163  INTEGER em_ifc(klev)
     164
     165  REAL em_pradj
     166  REAL em_cldf(klev), em_cldq(klev)
     167  REAL em_ftadj(klev), em_fradj(klev)
     168
     169  INTEGER ifc(klon, klev)
     170  REAL pradj(klon)
     171  REAL cldf(klon, klev), cldq(klon, klev)
     172  REAL ftadj(klon, klev), fqadj(klon, klev)
     173
     174  ! sb --
     175
     176  ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     177
     178  include "YOMCST.h"
     179  include "YOETHF.h"
     180  include "FCTTRE.h"
     181
     182  IF (first) THEN
     183
     184    ALLOCATE (em_t(klev))
     185    ALLOCATE (em_q(klev))
     186    ALLOCATE (em_qs(klev))
     187    ALLOCATE (em_u(klev), em_v(klev), em_tra(klev,nbtr))
     188    ALLOCATE (em_ph(klev+1), em_p(klev))
     189    ALLOCATE (em_work1(klev), em_work2(klev))
     190    ALLOCATE (em_d_t(klev), em_d_q(klev))
     191    ALLOCATE (em_d_u(klev), em_d_v(klev), em_d_tra(klev,nbtr))
     192    ALLOCATE (em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev))
     193    ALLOCATE (emmip(klev))
     194    ALLOCATE (emmke(klev))
     195    ALLOCATE (emma(klev))
     196
     197    first = .FALSE.
     198  END IF
     199
     200  qcond_incld(:, :) = 0.
     201
     202  ! @$$      print*,'debut conema'
     203
     204  DO i = 1, klon
     205    DO l = 1, klev + 1
     206      em_ph(l) = paprs(i, l)/100.0
     207    END DO
     208
     209    DO l = 1, klev
     210      em_p(l) = pplay(i, l)/100.0
     211      em_t(l) = t(i, l)
     212      em_q(l) = q(i, l)
     213      em_u(l) = u(i, l)
     214      em_v(l) = v(i, l)
     215      DO itra = 1, ntra
     216        em_tra(l, itra) = tra(i, l, itra)
     217      END DO
     218      ! @$$      print*,'em_t',em_t
     219      ! @$$      print*,'em_q',em_q
     220      ! @$$      print*,'em_qs',em_qs
     221      ! @$$      print*,'em_u',em_u
     222      ! @$$      print*,'em_v',em_v
     223      ! @$$      print*,'em_tra',em_tra
     224      ! @$$      print*,'em_p',em_p
     225
     226
     227
     228      zx_t = em_t(l)
     229      zdelta = max(0., sign(1.,rtt-zx_t))
     230      zx_qs = r2es*foeew(zx_t, zdelta)/em_p(l)/100.0
     231      zx_qs = min(0.5, zx_qs)
     232      ! @$$       print*,'zx_qs',zx_qs
     233      zcor = 1./(1.-retv*zx_qs)
     234      zx_qs = zx_qs*zcor
     235      em_qs(l) = zx_qs
     236      ! @$$      print*,'em_qs',em_qs
     237
     238      em_work1(l) = work1(i, l)
     239      em_work2(l) = work2(i, l)
     240      emmke(l) = 0
     241      ! emMa(l)=0
     242      ! Ma(i,l)=0
     243
     244      em_dtvpdt1(l) = 0.
     245      em_dtvpdq1(l) = 0.
     246      dtvpdt1(i, l) = 0.
     247      dtvpdq1(i, l) = 0.
     248    END DO
     249
     250    em_dplcldt = 0.
     251    em_dplcldr = 0.
     252    rain(i) = 0.0
     253    snow(i) = 0.0
     254    kbas(i) = 1
     255    ktop(i) = 1
     256    ! ajout SB:
     257    bas = 1
     258    top = 1
     259
     260
     261    ! sb3d      write(*,1792) (em_work1(m),m=1,klev)
     2621792 FORMAT ('sig avant convect ', /, 10(1X,E13.5))
     263
     264    ! sb d      write(*,1793) (em_work2(m),m=1,klev)
     2651793 FORMAT ('w avant convect ', /, 10(1X,E13.5))
     266
     267    ! @$$      print*,'avant convect'
     268    ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     269
     270
     271    ! print*,'avant convect i=',i
     272    CALL convect3(dtime, epmax, ok_adj_ema, em_t, em_q, em_qs, em_u, em_v, &
     273      em_tra, em_p, em_ph, klev, klev+1, klev-1, ntra, dtime, iflag, em_d_t, &
     274      em_d_q, em_d_u, em_d_v, em_d_tra, em_precip, em_bas, em_top, em_upwd, &
     275      em_dnwd, em_dnwdbis, em_work1, em_work2, emmip, emmke, emma, ment, &
     276      qent, tps, tls, sij, em_cape, em_tvp, em_pbase, em_bbase, em_dtvpdt1, &
     277      em_dtvpdq1, em_dplcldt, em_dplcldr, & ! sbl
     278      em_d_t2, em_d_q2, em_d_u2, em_d_v2, em_wd, em_qcond, em_qcondc) !sbl
     279    ! print*,'apres convect '
     280
     281    ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     282
     283    ! -- sb: Appel schema statistique de nuages couple a la convection
     284    ! (Bony et Emanuel 2001):
     285
     286    ! -- creer cvthermo.h qui contiendra les cstes thermo de LMDZ:
     287
     288    iflagcon = 3
     289    ! CALL cv_thermo(iflagcon)
     290
     291    ! -- appel schema de nuages:
     292
     293    ! CALL CLOUDS_SUB_LS(klev,em_q,em_qs,em_t
     294    ! i          ,em_p,em_ph,dtime,em_qcondc
     295    ! o          ,em_cldf,em_cldq,em_pradj,em_ftadj,em_fradj,em_ifc)
     296
     297    DO k = 1, klev
     298      cldf(i, k) = em_cldf(k) ! cloud fraction (0-1)
     299      cldq(i, k) = em_cldq(k) ! in-cloud water content (kg/kg)
     300      ftadj(i, k) = em_ftadj(k) ! (dT/dt)_{LS adj} (K/s)
     301      fqadj(i, k) = em_fradj(k) ! (dq/dt)_{LS adj} (kg/kg/s)
     302      ifc(i, k) = em_ifc(k) ! flag convergence clouds_gno (1 ou 2)
     303    END DO
     304    pradj(i) = em_pradj ! precip from LS supersat adj (mm/day)
     305
     306    ! sb --
     307
     308    ! SB:
     309    IF (iflag/=1 .AND. iflag/=4) THEN
     310      em_cape = 0.
    211311      DO l = 1, klev
    212          em_p(l) = pplay(i,l) / 100.0
    213          em_t(l) = t(i,l)
    214          em_q(l) = q(i,l)
    215          em_u(l) = u(i,l)
    216          em_v(l) = v(i,l)
    217          do itra = 1, ntra
    218           em_tra(l,itra) = tra(i,l,itra)
    219          enddo
    220 c@$$      print*,'em_t',em_t
    221 c@$$      print*,'em_q',em_q
    222 c@$$      print*,'em_qs',em_qs
    223 c@$$      print*,'em_u',em_u
    224 c@$$      print*,'em_v',em_v
    225 c@$$      print*,'em_tra',em_tra
    226 c@$$      print*,'em_p',em_p
    227 
    228  
    229 c
    230          zx_t = em_t(l)
    231          zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
    232          zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(l)/100.0
    233          zx_qs=MIN(0.5,zx_qs)
    234 c@$$       print*,'zx_qs',zx_qs
    235          zcor=1./(1.-retv*zx_qs)
    236          zx_qs=zx_qs*zcor
    237          em_qs(l) = zx_qs
    238 c@$$      print*,'em_qs',em_qs
    239 c
    240          em_work1(l) = work1(i,l)
    241          em_work2(l) = work2(i,l)
    242          emMke(l)=0
    243 c        emMa(l)=0
    244 c        Ma(i,l)=0
    245      
    246          em_dtvpdt1(l) = 0.
    247          em_dtvpdq1(l) = 0.
    248          dtvpdt1(i,l) = 0.
    249          dtvpdq1(i,l) = 0.
    250       ENDDO
    251 c
    252       em_dplcldt = 0.
    253       em_dplcldr = 0.
    254       rain(i) = 0.0
    255       snow(i) = 0.0
    256       kbas(i) = 1
    257       ktop(i) = 1
    258 c ajout SB:
    259       bas = 1
    260       top = 1
    261  
    262  
    263 c sb3d      write(*,1792) (em_work1(m),m=1,klev)
    264 1792  format('sig avant convect ',/,10(1X,E13.5))
    265 c
    266 c sb d      write(*,1793) (em_work2(m),m=1,klev)
    267 1793  format('w avant convect ',/,10(1X,E13.5))
    268  
    269 c@$$      print*,'avant convect'
    270 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    271 c
    272 
    273 c     print*,'avant convect i=',i
    274       CALL convect3(dtime,epmax,ok_adj_ema,
    275      .              em_t, em_q, em_qs,em_u ,em_v ,
    276      .              em_tra, em_p, em_ph,
    277      .              klev, klev+1, klev-1,ntra, dtime, iflag,
    278      .              em_d_t, em_d_q,em_d_u,em_d_v,
    279      .              em_d_tra, em_precip,
    280      .              em_bas, em_top,em_upwd, em_dnwd, em_dnwdbis,
    281      .              em_work1, em_work2,emmip,emMke,emMa,Ment,
    282      .  Qent,TPS,TLS,SIJ,em_CAPE,em_TVP,em_pbase,em_bbase,
    283      .  em_dtvpdt1,em_dtvpdq1,em_dplcldt,em_dplcldr, ! sbl
    284      .  em_d_t2,em_d_q2,em_d_u2,em_d_v2,em_wd,em_qcond,em_qcondc)!sbl
    285 c     print*,'apres convect '
    286 c
    287 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    288 c
    289 c -- sb: Appel schema statistique de nuages couple a la convection
    290 c (Bony et Emanuel 2001):
    291 
    292 c -- creer cvthermo.h qui contiendra les cstes thermo de LMDZ:
    293 
    294         iflagcon = 3
    295 c       CALL cv_thermo(iflagcon)
    296 
    297 c -- appel schema de nuages:
    298 
    299 c       CALL CLOUDS_SUB_LS(klev,em_q,em_qs,em_t
    300 c    i          ,em_p,em_ph,dtime,em_qcondc
    301 c    o          ,em_cldf,em_cldq,em_pradj,em_ftadj,em_fradj,em_ifc)
    302 
    303         do k = 1, klev
    304          cldf(i,k)  = em_cldf(k)  ! cloud fraction (0-1)
    305          cldq(i,k)  = em_cldq(k)  ! in-cloud water content (kg/kg)
    306          ftadj(i,k) = em_ftadj(k) ! (dT/dt)_{LS adj} (K/s)
    307          fqadj(i,k) = em_fradj(k) ! (dq/dt)_{LS adj} (kg/kg/s)
    308          ifc(i,k)   = em_ifc(k)   ! flag convergence clouds_gno (1 ou 2)
    309         enddo
    310         pradj(i) = em_pradj       ! precip from LS supersat adj (mm/day)
    311 
    312 c sb --
    313 c
    314 c SB:
    315       if (iflag.ne.1 .and. iflag.ne.4) then
    316          em_CAPE = 0.
    317       do l = 1, klev
    318          em_upwd(l) = 0.
    319          em_dnwd(l) = 0.
    320          em_dnwdbis(l) = 0.
    321          emMa(l) = 0.
    322          em_TVP(l) = 0.
    323       enddo
    324       endif
    325 c fin SB
    326 c
    327 c  If sig has been set to zero, then set Ma to zero
    328 c
    329       sigsum = 0.
    330       do k = 1,klev
    331         sigsum = sigsum + em_work1(k)
    332       enddo
    333       if (sigsum .eq. 0.0) then
    334         do k = 1,klev
    335           emMa(k) = 0.
    336         enddo
    337       endif
    338 c
    339 c sb3d       print*,'i, iflag=',i,iflag
    340 c
    341 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    342 c
    343 c       SORTIE DES ICB ET INB
    344 c       en fait inb et icb correspondent au niveau ou se trouve
    345 c       le nuage,le numero d'interface
    346 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    347  
    348 c modif SB:
    349       if (iflag.EQ.1 .or. iflag.EQ.4) then
    350        top=em_top
    351        bas=em_bas
    352        kbas(i) = em_bas
    353        ktop(i) = em_top
    354       endif
    355  
    356       pbase(i) = em_pbase
    357       bbase(i) = em_bbase
    358       rain(i) = em_precip/ 86400.0
    359       snow(i) = 0.0
    360       cape(i) = em_CAPE
    361       wd(i) = em_wd
    362       rflag(i) = REAL(iflag)
    363 c SB      kbas(i) = em_bas
    364 c SB      ktop(i) = em_top
    365       dplcldt(i) = em_dplcldt
    366       dplcldr(i) = em_dplcldr
    367       DO l = 1, klev
    368          d_t2(i,l) = dtime * em_d_t2(l)
    369          d_q2(i,l) = dtime * em_d_q2(l)
    370          d_u2(i,l) = dtime * em_d_u2(l)
    371          d_v2(i,l) = dtime * em_d_v2(l)
    372 
    373          d_t(i,l) = dtime * em_d_t(l)
    374          d_q(i,l) = dtime * em_d_q(l)
    375          d_u(i,l) = dtime * em_d_u(l)
    376          d_v(i,l) = dtime * em_d_v(l)
    377          do itra = 1, ntra
    378          d_tra(i,l,itra) = dtime * em_d_tra(l,itra)
    379          enddo
    380          upwd(i,l) = em_upwd(l)
    381          dnwd(i,l) = em_dnwd(l)
    382          dnwdbis(i,l) = em_dnwdbis(l)
    383          work1(i,l) = em_work1(l)
    384          work2(i,l) = em_work2(l)
    385          Ma(i,l)=emMa(l)
    386          tvp(i,l)=em_TVP(l)
    387          dtvpdt1(i,l) = em_dtvpdt1(l)
    388          dtvpdq1(i,l) = em_dtvpdq1(l)
    389 
    390          if (iflag_clw.eq.0) then
    391             qcond_incld(i,l) = em_qcondc(l)
    392          else if (iflag_clw.eq.1) then
    393             qcond_incld(i,l) = em_qcond(l)
    394          endif
    395       ENDDO
    396   999 CONTINUE
    397 
    398 c   On calcule une eau liquide diagnostique en fonction de la
    399 c  precip.
    400       if ( iflag_clw.eq.2 ) then
    401       do l=1,klev
    402          do i=1,klon
    403             if (ktop(i)-kbas(i).gt.0.and.
    404      s         l.ge.kbas(i).and.l.le.ktop(i)) then
    405                qcond_incld(i,l)=rain(i)*8.e4
    406 c    s         *(pplay(i,l      )-paprs(i,ktop(i)+1))
    407      s         /(pplay(i,kbas(i))-pplay(i,ktop(i)))
    408 c    s         **2
    409             else
    410                qcond_incld(i,l)=0.
    411             endif
    412          enddo
    413          print*,'l=',l,',   qcond_incld=',qcond_incld(1,l)
    414       enddo
    415       endif
    416  
    417 
    418       RETURN
    419       END
    420 
     312        em_upwd(l) = 0.
     313        em_dnwd(l) = 0.
     314        em_dnwdbis(l) = 0.
     315        emma(l) = 0.
     316        em_tvp(l) = 0.
     317      END DO
     318    END IF
     319    ! fin SB
     320
     321    ! If sig has been set to zero, then set Ma to zero
     322
     323    sigsum = 0.
     324    DO k = 1, klev
     325      sigsum = sigsum + em_work1(k)
     326    END DO
     327    IF (sigsum==0.0) THEN
     328      DO k = 1, klev
     329        emma(k) = 0.
     330      END DO
     331    END IF
     332
     333    ! sb3d       print*,'i, iflag=',i,iflag
     334
     335    ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     336
     337    ! SORTIE DES ICB ET INB
     338    ! en fait inb et icb correspondent au niveau ou se trouve
     339    ! le nuage,le numero d'interface
     340    ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     341
     342    ! modif SB:
     343    IF (iflag==1 .OR. iflag==4) THEN
     344      top = em_top
     345      bas = em_bas
     346      kbas(i) = em_bas
     347      ktop(i) = em_top
     348    END IF
     349
     350    pbase(i) = em_pbase
     351    bbase(i) = em_bbase
     352    rain(i) = em_precip/86400.0
     353    snow(i) = 0.0
     354    cape(i) = em_cape
     355    wd(i) = em_wd
     356    rflag(i) = real(iflag)
     357    ! SB      kbas(i) = em_bas
     358    ! SB      ktop(i) = em_top
     359    dplcldt(i) = em_dplcldt
     360    dplcldr(i) = em_dplcldr
     361    DO l = 1, klev
     362      d_t2(i, l) = dtime*em_d_t2(l)
     363      d_q2(i, l) = dtime*em_d_q2(l)
     364      d_u2(i, l) = dtime*em_d_u2(l)
     365      d_v2(i, l) = dtime*em_d_v2(l)
     366
     367      d_t(i, l) = dtime*em_d_t(l)
     368      d_q(i, l) = dtime*em_d_q(l)
     369      d_u(i, l) = dtime*em_d_u(l)
     370      d_v(i, l) = dtime*em_d_v(l)
     371      DO itra = 1, ntra
     372        d_tra(i, l, itra) = dtime*em_d_tra(l, itra)
     373      END DO
     374      upwd(i, l) = em_upwd(l)
     375      dnwd(i, l) = em_dnwd(l)
     376      dnwdbis(i, l) = em_dnwdbis(l)
     377      work1(i, l) = em_work1(l)
     378      work2(i, l) = em_work2(l)
     379      ma(i, l) = emma(l)
     380      tvp(i, l) = em_tvp(l)
     381      dtvpdt1(i, l) = em_dtvpdt1(l)
     382      dtvpdq1(i, l) = em_dtvpdq1(l)
     383
     384      IF (iflag_clw==0) THEN
     385        qcond_incld(i, l) = em_qcondc(l)
     386      ELSE IF (iflag_clw==1) THEN
     387        qcond_incld(i, l) = em_qcond(l)
     388      END IF
     389    END DO
     390  END DO
     391
     392  ! On calcule une eau liquide diagnostique en fonction de la
     393  ! precip.
     394  IF (iflag_clw==2) THEN
     395    DO l = 1, klev
     396      DO i = 1, klon
     397        IF (ktop(i)-kbas(i)>0 .AND. l>=kbas(i) .AND. l<=ktop(i)) THEN
     398          qcond_incld(i, l) = rain(i)*8.E4 & ! s         *(pplay(i,l
     399                                             ! )-paprs(i,ktop(i)+1))
     400            /(pplay(i,kbas(i))-pplay(i,ktop(i)))
     401          ! s         **2
     402        ELSE
     403          qcond_incld(i, l) = 0.
     404        END IF
     405      END DO
     406      PRINT *, 'l=', l, ',   qcond_incld=', qcond_incld(1, l)
     407    END DO
     408  END IF
     409
     410
     411  RETURN
     412END SUBROUTINE conema3
     413
Note: See TracChangeset for help on using the changeset viewer.