Changeset 4654 for LMDZ6


Ignore:
Timestamp:
Aug 30, 2023, 9:30:12 AM (9 months ago)
Author:
fhourdin
Message:

Replayisation de lscp_mod et vdif_kcay.F90

Location:
LMDZ6/trunk/libf/phylmd
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90

    r4230 r4654  
    167167!   iflag_pbl peut etre utilise comme longuer de melange
    168168       IF (iflag_pbl.GE.31) THEN
    169           CALL vdif_kcay(knon,dtime,RG,RD,ypaprs,yt, &
     169          if (knon>0) then
     170          CALL vdif_kcay(knon,klev,dtime,RG,RD,ypaprs,yt, &
    170171               yzlev,yzlay,yu,yv,yteta, &
    171172               ycdragm,yq2,q2diag,ykmm,ykmn,yustar, &
    172173               iflag_pbl)
     174          endif
    173175       ELSE IF (iflag_pbl<20) THEN
    174176          CALL yamada4(ni,nsrf,knon,dtime,RG,RD,ypaprs,yt, &

  • LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90

    r4651 r4654  
    1 module lscp_ini_mod
    2 
    3 implicit none
     1MODULE lscp_ini_mod
     2
     3IMPLICIT NONE
    44
    55  ! PARAMETERS for lscp:
    … …  
    259259    CALL cloudth_ini(iflag_cloudth_vert,iflag_ratqs)
    260260
    261 end subroutine lscp_ini
    262 
    263 end module lscp_ini_mod
     261RETURN
     262
     263END SUBROUTINE lscp_ini
     264
     265END MODULE lscp_ini_mod

  • LMDZ6/trunk/libf/phylmd/lscp_mod.F90

    r4651 r4654  
    77!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    88SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
    9      paprs,pplay,t,q,ptconv,ratqs,                      &
     9     paprs,pplay,temp,q,ptconv,ratqs,                      &
    1010     d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri,  &
    1111     pfraclr,pfracld,                                   &
    … …  
    9191!
    9292USE print_control_mod, ONLY: prt_level, lunout
     93
     94! USE de modules contenant des fonctions.
    9395USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
    9496USE lscp_tools_mod, ONLY : calc_qsat_ecmwf, icefrac_lscp, calc_gammasat
    … …  
    9698USE ice_sursat_mod, ONLY : ice_sursat
    9799
     100! Use du module lscp_ini_mod contenant les constantes
    98101USE lscp_ini_mod, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
    99102USE lscp_ini_mod, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
    … …  
    119122  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
    120123  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
    121   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: t               ! temperature (K)
     124  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
    122125  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: q               ! total specific humidity (= vapor phase) [kg/kg]
    123126  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
    … …  
    194197  !----------------
    195198
    196   REAL qsl(klon), qsi(klon)
    197   REAL zct, zcl,zexpo
    198   REAL ctot(klon,klev)
    199   REAL ctot_vol(klon,klev)
    200   REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
    201   REAL zdqsdT_raw(klon)
    202   REAL gammasat(klon),dgammasatdt(klon)                ! coefficient to make cold condensation at the correct RH and derivative wrt T
    203   REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
    204   REAL cste
    205   REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
    206   REAL Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
    207   REAL erf
    208   REAL qcloud(klon), icefrac_mpc(klon,klev), qincloud_mpc(klon)
    209   REAL zrfl(klon), zrfln(klon), zqev, zqevt
    210   REAL zifl(klon), zifln(klon), ziflprev(klon),zqev0,zqevi, zqevti
    211   REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon)
    212   REAL zoliql(klon), zoliqi(klon)
    213   REAL zt(klon)
    214   REAL zdz(klon),zrho(klon,klev),iwc(klon)
    215   REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon)
    216   REAL zmelt,zrain,zsnow,zprecip
    217   REAL dzfice(klon)
    218   REAL zsolid
    219   REAL qtot(klon), qzero(klon)
    220   REAL dqsl(klon),dqsi(klon)
    221   REAL smallestreal
     199  REAL,DIMENSION(klon) :: qsl, qsi
     200  REAL :: zct, zcl,zexpo
     201  REAL, DIMENSION(klon,klev) :: ctot
     202  REAL, DIMENSION(klon,klev) :: ctot_vol
     203  REAL, DIMENSION(klon) :: zqs, zdqs
     204  REAL :: zdelta, zcor, zcvm5
     205  REAL, DIMENSION(klon) :: zdqsdT_raw
     206  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                ! coefficient to make cold condensation at the correct RH and derivative wrt T
     207  REAL, DIMENSION(klon) :: Tbef,qlbef,DT
     208  REAL :: num,denom
     209  REAL :: cste
     210  REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta
     211  REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2
     212  REAL :: erf
     213  REAL, DIMENSION(klon,klev) :: icefrac_mpc
     214  REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
     215  REAL, DIMENSION(klon) :: zrfl, zrfln
     216  REAL :: zqev, zqevt
     217  REAL, DIMENSION(klon) :: zifl, zifln, ziflprev
     218  REAL :: zqev0,zqevi, zqevti
     219  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
     220  REAL, DIMENSION(klon) :: zoliql, zoliqi
     221  REAL, DIMENSION(klon) :: zt
     222  REAL, DIMENSION(klon,klev) :: zrho
     223  REAL, DIMENSION(klon) :: zdz,iwc
     224  REAL :: zchau,zfroi
     225  REAL, DIMENSION(klon) :: zfice,zneb,znebprecip
     226  REAL :: zmelt,zrain,zsnow,zprecip
     227  REAL, DIMENSION(klon) :: dzfice
     228  REAL :: zsolid
     229  REAL, DIMENSION(klon) :: qtot, qzero
     230  REAL, DIMENSION(klon) :: dqsl,dqsi
     231  REAL :: smallestreal
    222232  !  Variables for Bergeron process
    223   REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
    224   REAL zqpreci(klon), zqprecl(klon)
     233  REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
     234  REAL, DIMENSION(klon) :: zqpreci, zqprecl
    225235  ! Variables precipitation energy conservation
    226   REAL zmqc(klon)
    227   REAL zalpha_tr
    228   REAL zfrac_lessi
    229   REAL zprec_cond(klon)
    230   REAL zmair(klon), zcpair, zcpeau
    231   REAL zlh_solid(klon), zm_solid         ! for liquid -> solid conversion
    232   REAL zrflclr(klon), zrflcld(klon)
    233   REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon)
    234   REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon)
    235   REAL ziflclr(klon), ziflcld(klon)
    236   REAL znebprecipclr(klon), znebprecipcld(klon)
    237   REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon)
    238   REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon)
    239   REAL velo(klon,klev), vr, ffallv
    240   REAL qlmpc, qimpc, rnebmpc
    241   REAL radocondi(klon,klev), radocondl(klon,klev)
    242   REAL effective_zneb
    243   REAL distcltop1D(klon), temp_cltop1D(klon)
    244 
    245 
    246   INTEGER i, k, n, kk
    247   INTEGER n_i(klon), iter
     236  REAL, DIMENSION(klon) :: zmqc
     237  REAL :: zalpha_tr
     238  REAL :: zfrac_lessi
     239  REAL, DIMENSION(klon) :: zprec_cond
     240  REAL, DIMENSION(klon) :: zmair
     241  REAL :: zcpair, zcpeau
     242  REAL, DIMENSION(klon) :: zlh_solid
     243  REAL :: zm_solid         ! for liquid -> solid conversion
     244  REAL, DIMENSION(klon) :: zrflclr, zrflcld
     245  REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
     246  REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
     247  REAL, DIMENSION(klon) :: ziflclr, ziflcld
     248  REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld
     249  REAL, DIMENSION(klon) :: tot_zneb, tot_znebn, d_tot_zneb
     250  REAL, DIMENSION(klon) :: d_znebprecip_clr_cld, d_znebprecip_cld_clr
     251  REAL, DIMENSION(klon,klev) :: velo
     252  REAL :: vr, ffallv
     253  REAL :: qlmpc, qimpc, rnebmpc
     254  REAL, DIMENSION(klon,klev) :: radocondi, radocondl
     255  REAL :: effective_zneb
     256  REAL, DIMENSION(klon) :: distcltop1D, temp_cltop1D
     257
     258
     259  INTEGER i, k, n, kk, iter
     260  INTEGER, DIMENSION(klon) :: n_i
    248261  INTEGER ncoreczq
    249   INTEGER mpc_bl_points(klon,klev)
    250 
    251   LOGICAL lognormale(klon)
    252   LOGICAL keepgoing(klon)
    253   LOGICAL,SAVE :: appel1er
    254   !$OMP THREADPRIVATE(appel1er)
     262  INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
     263
     264  LOGICAL, DIMENSION(klon) :: lognormale
     265  LOGICAL, DIMENSION(klon) :: keepgoing
    255266
    256267  CHARACTER (len = 20) :: modname = 'lscp'
    257268  CHARACTER (len = 80) :: abort_message
    258269 
    259   DATA appel1er /.TRUE./
    260270
    261271!===============================================================================
    … …  
    264274
    265275! Few initial checks
     276
    266277
    267278IF (iflag_fisrtilp_qsat .LT. 0) THEN
    … …  
    332343fcontrN(:,:) = 0.0
    333344fcontrP(:,:) = 0.0
     345distcltop(:,:)=0.
     346temp_cltop(:,:)=0.
    334347
    335348!c_iso: variable initialisation for iso
    … …  
    346359    ! Initialisation temperature and specific humidity
    347360    DO i = 1, klon
    348         zt(i)=t(i,k)
     361        zt(i)=temp(i,k)
    349362        zq(i)=q(i,k)
    350363        !c_iso init of iso
    … …  
    382395            zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
    383396            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
    384             zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
     397            zt(i) = ( (temp(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
    385398             / (zcpair + zmqc(i)*zcpeau)
    386399
    … …  
    625638                         zq,zqta,fraca,                            &
    626639                         qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
    627                          ratqs,zqs,t,                              &
     640                         ratqs,zqs,temp,                              &
    628641                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
    629642
    … …  
    637650                         zq,zqta,fraca,                                     &
    638651                         qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
    639                          ratqs,zqs,t, &
     652                         ratqs,zqs,temp, &
    640653                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
    641654
    … …  
    648661                         zq,zqta,fraca,                                     &
    649662                         qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
    650                          ratqs,zqs,t, &
     663                         ratqs,zqs,temp, &
    651664                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
    652665
    … …  
    658671
    659672                    CALL cloudth_mpc(klon,klev,k,iflag_mpc_bl,mpc_bl_points, &
    660                          t,ztv,zq,zqta,fraca, zpspsk,                      &
     673                         temp,ztv,zq,zqta,fraca, zpspsk,                      &
    661674                         paprs,pplay,ztla,zthl,ratqs,zqs,psfl,             &
    662675                         qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol,    &
    … …  
    739752                  IF (iflag_t_glace.GE.4) THEN
    740753                  ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
    741                        CALL distance_to_cloud_top(klon,klev,k,t,pplay,paprs,rneb,distcltop1D,temp_cltop1D)
     754                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D)
    742755                  ENDIF
    743756                  CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, distcltop1D(:),temp_cltop1D(:),zfice(:),dzfice(:))
    … …  
    780793                        !------------------------------------
    781794
    782                         CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), &
     795                        CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, temp(i,k), zq(i), &
    783796                             gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),                 &
    784797                             rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                   &
    … …  
    837850        ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
    838851        IF (iflag_t_glace.GE.4) THEN
    839             CALL distance_to_cloud_top(klon,klev,k,t,pplay,paprs,rneb,distcltop1D,temp_cltop1D)
     852            CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D)
    840853            distcltop(:,k)=distcltop1D(:)
    841854            temp_cltop(:,k)=temp_cltop1D(:)
    … …  
    11961209        d_q(i,k) = zq(i) - q(i,k)
    11971210        ! c_iso: same for isotopes
    1198         d_t(i,k) = zt(i) - t(i,k)
     1211        d_t(i,k) = zt(i) - temp(i,k)
    11991212    ENDDO
    12001213
    … …  
    12451258        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
    12461259
    1247             IF (t(i,k) .GE. t_glace_min) THEN
     1260            IF (temp(i,k) .GE. t_glace_min) THEN
    12481261                zalpha_tr = a_tr_sca(3)
    12491262            ELSE
    … …  
    12691282            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
    12701283
    1271                 IF (t(i,kk) .GE. t_glace_min) THEN
     1284                IF (temp(i,kk) .GE. t_glace_min) THEN
    12721285                    zalpha_tr = a_tr_sca(1)
    12731286                ELSE
    … …  
    13181331  ENDIF
    13191332
     1333
     1334RETURN
     1335
    13201336END SUBROUTINE lscp
    13211337!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  • LMDZ6/trunk/libf/phylmd/physiqex_mod.F90

    r4652 r4654  
    99     &            debut,lafin,pdtphys, &
    1010     &            paprs,pplay,pphi,pphis,presnivs, &
    11      &            u,v,rot,t,qx, &
     11     &            u,v,rot,temp,qx, &
    1212     &            flxmass_w, &
    1313     &            d_u, d_v, d_t, d_qx, d_ps)
    1414
    15       USE dimphy, only : klon,klev
    16       USE infotrac_phy, only : nqtot
    17       USE geometry_mod, only : latitude
    18 !      USE comcstphy, only : rg
    19       USE ioipsl, only : ymds2ju
    20       USE phys_state_var_mod, only : phys_state_var_init
    21       USE phyetat0_mod, only: phyetat0
    22       USE output_physiqex_mod, ONLY: output_physiqex
     15
     16USE dimphy, only : klon,klev
     17USE infotrac_phy, only : nqtot
     18USE geometry_mod, only : latitude
     19USE ioipsl, only : ymds2ju
     20USE phys_state_var_mod, only : phys_state_var_init
     21USE phyetat0_mod, only: phyetat0
     22USE output_physiqex_mod, ONLY: output_physiqex
     23use vdif_ini, only : vdif_ini_
     24USE lmdz_thermcell_ini, ONLY : thermcell_ini
     25USE ioipsl_getin_p_mod, ONLY : getin_p
     26USE wxios, ONLY: missing_val, using_xios
     27USE lscp_mod, ONLY : lscp
     28USE lscp_ini_mod, ONLY : lscp_ini
     29USE add_phys_tend_mod, ONLY : add_phys_tend
     30
    2331
    2432      IMPLICIT none
     33
     34include "YOETHF.h"
     35
     36
     37
     38
    2539!
    2640! Routine argument:
    … …  
    3246      logical,intent(in) :: lafin ! signals last call to physics
    3347      real,intent(in) :: pdtphys ! physics time step (s)
    34       real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa)
    35       real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa)
    36       real,intent(in) :: pphi(klon,klev) ! geopotential at mid-layer
    37       real,intent(in) :: pphis(klon) ! surface geopotential
    38       real,intent(in) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers
    39       real,intent(in) :: u(klon,klev) ! eastward zonal wind (m/s)
    40       real,intent(in) :: v(klon,klev) ! northward meridional wind (m/s)
    41       real,intent(in) :: rot(klon,klev) ! northward meridional wind (m/s)
    42       real,intent(in) :: t(klon,klev) ! temperature (K)
    43       real,intent(in) :: qx(klon,klev,nqtot) ! tracers (.../kg_air)
    44       real,intent(in) :: flxmass_w(klon,klev) ! vertical mass flux
    45       real,intent(out) :: d_u(klon,klev) ! physics tendency on u (m/s/s)
    46       real,intent(out) :: d_v(klon,klev) ! physics tendency on v (m/s/s)
    47       real,intent(out) :: d_t(klon,klev) ! physics tendency on t (K/s)
    48       real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers
    49       real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure
     48      real,dimension(klon,klev+1),intent(in) :: paprs ! interlayer pressure (Pa)
     49      real,dimension(klon,klev),intent(in) :: pplay ! mid-layer pressure (Pa)
     50      real,dimension(klon,klev),intent(in) :: pphi ! geopotential at mid-layer
     51      real,dimension(klon),intent(in) :: pphis ! surface geopotential
     52      real,dimension(klev),intent(in) :: presnivs ! pseudo-pressure (Pa) of mid-layers
     53      real,dimension(klon,klev),intent(in) :: u ! eastward zonal wind (m/s)
     54      real,dimension(klon,klev),intent(in) :: v ! northward meridional wind (m/s)
     55      real,dimension(klon,klev),intent(in) :: rot ! northward meridional wind (m/s)
     56      real,dimension(klon,klev),intent(in) :: temp ! temperature (K)
     57      real,dimension(klon,klev,nqtot),intent(in) :: qx  ! tracers (.../kg_air)
     58      real,dimension(klon,klev),intent(in) :: flxmass_w ! vertical mass flux
     59      real,dimension(klon,klev),intent(out) :: d_u ! physics tendency on u (m/s/s)
     60      real,dimension(klon,klev),intent(out) :: d_v ! physics tendency on v (m/s/s)
     61      real,dimension(klon,klev),intent(out) :: d_t ! physics tendency on t (K/s)
     62      real,dimension(klon,klev,nqtot),intent(out) :: d_qx ! physics tendency on tracers
     63      real,dimension(klon),intent(out) :: d_ps ! physics tendency on surface pressure
     64
     65      real, dimension(klon,klev) :: u_loc
     66      real, dimension(klon,klev) :: v_loc
     67      real, dimension(klon,klev) :: t_loc
     68      real, dimension(klon,klev) :: h_loc
     69      real, dimension(klon,klev) :: d_u_loc,d_v_loc,d_t_loc,d_h_loc
     70
     71      real, dimension(klon,klev) :: d_u_dyn,d_v_dyn,d_t_dyn
     72      real, dimension(klon,klev,nqtot) :: d_q_dyn
     73      real, allocatable, dimension(:,:), save :: u_prev,v_prev,t_prev
     74      real, allocatable, dimension(:,:,:), save :: q_prev
     75!$OMP THREADPRIVATE(u_prev,v_prev,t_prev,q_prev)
     76
     77
     78
     79      real, dimension(klon,klev) :: d_u_vdif,d_v_vdif,d_t_vdif,d_h_vdif
     80      real, dimension(klon,klev) :: d_u_the,d_v_the,d_t_the
     81      real, dimension(klon,klev,nqtot) :: q_loc,d_q_loc,d_q_vdif,d_q_the
     82
     83real, dimension(klon) :: capcal,z0m,z0h,dtsrf,emis,fluxsrf,cdh,cdv,tsrf_
     84real, dimension(klon,klev) :: zzlay,masse,zpopsk
     85real, dimension(klon,klev+1) :: zzlev,kz_v,kz_h,richardson
     86
     87real, save, allocatable, dimension(:) :: tsrf,f0,zmax0
     88real, save, allocatable, dimension(:,:) :: q2
     89!$OMP THREADPRIVATE(tsrf,q2,f0,zmax0)
     90
     91    real,save :: ratqsbas=0.002,ratqshaut=0.3,ratqsp0=50000.,ratqsdp=20000.
     92    !$OMP THREADPRIVATE(ratqsbas,ratqshaut,ratqsp0,ratqsdp)
     93
     94
     95real :: z1,z2,tau_thermals
     96logical :: lwrite
     97integer :: iflag_replay
     98
     99integer :: iflag_thermals=18
     100
     101!--------------------------------------------------------------
     102! Declaration lscp
     103!--------------------------------------------------------------
     104  INTEGER                         :: iflag_cld_th    ! flag that determines the distribution of convective clouds    ! IN
     105  INTEGER                         :: iflag_ice_thermo! flag to activate the ice thermodynamics                       ! IN
     106  LOGICAL                         :: ok_ice_sursat   ! flag to determine if ice sursaturation is activated           ! IN
     107  LOGICAL, DIMENSION(klon,klev)   :: ptconv          ! grid points where deep convection scheme is active            ! IN
     108  REAL, DIMENSION(klon,klev)      :: ztv             ! virtual potential temperature [K]                             ! IN
     109  REAL, DIMENSION(klon,klev)      :: zqta            ! specific humidity within thermals [kg/kg]                     ! IN
     110  REAL, DIMENSION(klon,klev+1)      :: frac_the,fm_the
     111  REAL, DIMENSION(klon,klev)      :: zpspsk          ! exner potential (p/100000)**(R/cp)                            ! IN
     112  REAL, DIMENSION(klon,klev)      :: ztla            ! liquid temperature within thermals [K]                        ! IN
     113  REAL, DIMENSION(klon,klev)         :: zthl         ! liquid potential temperature [K]                              ! INOUT
     114  REAL, DIMENSION(klon,klev)      :: ratqs            ! function of pressure that sets the large-scale               ! INOUT
     115  REAL, DIMENSION(klon,klev)      :: beta             ! conversion rate of condensed water                           ! INOUT
     116  REAL, DIMENSION(klon,klev)      :: rneb_seri        ! fraction nuageuse en memoire                                 ! INOUT
     117  REAL, DIMENSION(klon,klev)      :: d_t_lscp         ! temperature increment [K]                                    ! OUT
     118  REAL, DIMENSION(klon,klev)      :: d_q_lscp         ! specific humidity increment [kg/kg]                          ! OUT
     119  REAL, DIMENSION(klon,klev)      :: d_ql_lscp        ! liquid water increment [kg/kg]                               ! OUT
     120  REAL, DIMENSION(klon,klev)      :: d_qi_lscp        ! cloud ice mass increment [kg/kg]                             ! OUT
     121  REAL, DIMENSION(klon,klev)      :: rneb             ! cloud fraction [-]                                           ! OUT
     122  REAL, DIMENSION(klon,klev)      :: rneblsvol        ! cloud fraction per unit volume [-]                           ! OUT
     123  REAL, DIMENSION(klon,klev)      :: pfraclr          ! precip fraction clear-sky part [-]                           ! OUT
     124  REAL, DIMENSION(klon,klev)      :: pfracld          ! precip fraction cloudy part [-]                              ! OUT
     125  REAL, DIMENSION(klon,klev)      :: radocond         ! condensed water used in the radiation scheme [kg/kg]         ! OUT
     126  REAL, DIMENSION(klon,klev)      :: radicefrac       ! ice fraction of condensed water for radiation scheme         ! OUT
     127  REAL, DIMENSION(klon,klev)      :: rhcl             ! clear-sky relative humidity [-]                              ! OUT
     128  REAL, DIMENSION(klon)           :: rain             ! surface large-scale rainfall [kg/s/m2]                       ! OUT
     129  REAL, DIMENSION(klon)           :: snow             ! surface large-scale snowfall [kg/s/m2]                       ! OUT
     130  REAL, DIMENSION(klon,klev)      :: qsatl            ! saturation specific humidity wrt liquid [kg/kg]              ! OUT
     131  REAL, DIMENSION(klon,klev)      :: qsats            ! saturation specific humidity wrt ice [kg/kg]                 ! OUT
     132  REAL, DIMENSION(klon,klev+1)    :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]            ! OUT
     133  REAL, DIMENSION(klon,klev+1)    :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]            ! OUT
     134  REAL, DIMENSION(klon,klev)      :: distcltop        ! distance to cloud top [m]                                    ! OUT
     135  REAL, DIMENSION(klon,klev)      :: temp_cltop       ! temperature of cloud top [K]                                 ! OUT
     136  REAL, DIMENSION(klon,klev)      :: frac_impa        ! scavenging fraction due tu impaction [-]                     ! OUT
     137  REAL, DIMENSION(klon,klev)      :: frac_nucl        ! scavenging fraction due tu nucleation [-]                    ! OUT
     138  REAL, DIMENSION(klon,klev)      :: qclr             ! specific total water content in clear sky region [kg/kg]     ! OUT
     139  REAL, DIMENSION(klon,klev)      :: qcld             ! specific total water content in cloudy region [kg/kg]        ! OUT
     140  REAL, DIMENSION(klon,klev)      :: qss              ! specific total water content in supersat region [kg/kg]      ! OUT
     141  REAL, DIMENSION(klon,klev)      :: qvc              ! specific vapor content in clouds [kg/kg]                     ! OUT
     142  REAL, DIMENSION(klon,klev)      :: rnebclr          ! mesh fraction of clear sky [-]                               ! OUT
     143  REAL, DIMENSION(klon,klev)      :: rnebss           ! mesh fraction of ISSR [-]                                    ! OUT
     144  REAL, DIMENSION(klon,klev)      :: gamma_ss         ! coefficient governing the ice nucleation RHi threshold [-]   ! OUT   
     145  REAL, DIMENSION(klon,klev)      :: Tcontr           ! threshold temperature for contrail formation [K]             ! OUT
     146  REAL, DIMENSION(klon,klev)      :: qcontr           ! threshold humidity for contrail formation [kg/kg]            ! OUT
     147  REAL, DIMENSION(klon,klev)      :: qcontr2          ! // (2nd expression more consistent with LMDZ expression of q)! OUT         
     148  REAL, DIMENSION(klon,klev)      :: fcontrN          ! fraction of grid favourable to non-persistent contrails      ! OUT
     149  REAL, DIMENSION(klon,klev)      :: fcontrP          ! fraction of grid favourable to persistent contrails          ! OUT
     150!--------------------------------------------------------------
     151
     152  REAL, DIMENSION(klon,klev)      :: d_t_eva,d_q_eva,d_ql_eva,d_qi_eva
     153include "YOMCST.h"
    50154
    51155!    include "clesphys.h"
    … …  
    58162
    59163
    60 real :: temp_newton(klon,klev)
    61 integer :: k
     164real,dimension(klon,klev) :: temp_newton
     165integer :: i,k,iq
     166INTEGER, SAVE :: itap=0
     167!$OMP THREADPRIVATE(itap)
     168INTEGER, SAVE :: abortphy=0   ! Reprere si on doit arreter en fin de phys
     169!$OMP THREADPRIVATE(abortphy)
     170
     171integer, save :: iflag_reevap=1,iflag_newton=0,iflag_vdif=1,iflag_lscp=1,iflag_cloudth_vert=3,iflag_ratqs=4
     172!$OMP THREADPRIVATE(iflag_reevap,iflag_newton,iflag_vdif,iflag_lscp,iflag_cloudth_vert,iflag_ratqs)
     173
    62174logical, save :: first=.true.
    63175!$OMP THREADPRIVATE(first)
    64 
    65 real,save :: rg=9.81
    66 !$OMP THREADPRIVATE(rg)
    67176
    68177! For I/Os
    69178integer :: itau0
    70179real :: zjulian
     180real,dimension(klon,klev) :: du0,dv0,dqbs0
     181real,dimension(klon,klev) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
    71182
    72183
    … …  
    93204  ! Initialize IOIPSL output file
    94205#endif
     206call suphel
     207call vdif_ini_(klon,RCPD,RD,RG,RKAPPA)
     208! Pourquoi ce tau_thermals en argument ??? AFAIRE
     209tau_thermals=0.
     210call getin_p('iflag_thermals',iflag_thermals)
     211
     212call getin_p('iflag_newton',iflag_newton)
     213call getin_p('iflag_reevap',iflag_reevap)
     214call getin_p('iflag_cloudth_vert',iflag_cloudth_vert)
     215call getin_p('iflag_ratqs',iflag_ratqs)
     216call getin_p('iflag_vdif',iflag_vdif)
     217call getin_p('iflag_lscp',iflag_lscp)
     218call getin_p('ratqsbas',ratqsbas)
     219call getin_p('ratqshaut',ratqshaut)
     220call getin_p('ratqsp0',ratqsp0)
     221call getin_p('ratqsdp',ratqsdp)
     222CALL thermcell_ini(iflag_thermals,0,tau_thermals,6, &
     223   &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
     224CALL lscp_ini(pdtphys,.false.,iflag_ratqs, RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG)
     225
     226
     227
     228allocate(tsrf(klon),q2(klon,klev+1),f0(klon),zmax0(klon))
     229allocate(u_prev(klon,klev),v_prev(klon,klev),t_prev(klon,klev),q_prev(klon,klev,nqtot))
     230
     231u_prev(:,:)=u(:,:)
     232v_prev(:,:)=v(:,:)
     233t_prev(:,:)=temp(:,:)
     234q_prev(:,:,:)=qx(:,:,:)
     235
     236q2=1.e-10
     237tsrf=temp(:,1)
     238f0=0.
     239zmax0=0.
     240
     241iflag_replay=0
     242call getin_p('iflag_replay',iflag_replay)
     243if ( iflag_replay >= 0 ) CALL iophys_ini(pdtphys)
     244
    95245
    96246endif ! of if (debut)
    … …  
    100250!------------------------------------------------------------
    101251
     252d_u_dyn(:,:)=(u(:,:)-u_prev(:,:))/pdtphys
     253d_v_dyn(:,:)=(v(:,:)-v_prev(:,:))/pdtphys
     254d_t_dyn(:,:)=(temp(:,:)-t_prev(:,:))/pdtphys
     255d_q_dyn(:,:,:)=(qx(:,:,:)-q_prev(:,:,:))/pdtphys
    102256
    103257! set all tendencies to zero
    … …  
    108262d_ps(1:klon)=0.
    109263
     264u_loc(1:klon,1:klev)=u(1:klon,1:klev)
     265v_loc(1:klon,1:klev)=v(1:klon,1:klev)
     266t_loc(1:klon,1:klev)=temp(1:klon,1:klev)
     267d_u_loc(1:klon,1:klev)=0.
     268d_v_loc(1:klon,1:klev)=0.
     269d_t_loc(1:klon,1:klev)=0.
     270do iq=1,nqtot
     271   do k=1,klev
     272      do i=1,klon
     273         q_loc(i,k,iq)=qx(i,k,iq)
     274      enddo
     275   enddo
     276enddo
     277
     278du0(1:klon,1:klev)=0.
     279dv0(1:klon,1:klev)=0.
     280dqbs0(1:klon,1:klev)=0.
     281
     282
     283
    110284!------------------------------------------------------------
    111285! Calculs
    112286!------------------------------------------------------------
    113287
    114 ! compute tendencies to return to the dynamics:
    115 ! "friction" on the first layer
    116 d_u(1:klon,1)=-u(1:klon,1)/86400.
    117 d_v(1:klon,1)=-v(1:klon,1)/86400.
    118 ! newtonian relaxation towards temp_newton()
     288!------------------------------------------------------------
     289! Rappel en temperature et frottement dans la premiere chouche
     290!------------------------------------------------------------
     291
     292if ( iflag_newton == 1 ) then
     293    ! compute tendencies to return to the dynamics:
     294    ! "friction" on the first layer
     295    d_u(1:klon,1)=-u(1:klon,1)/86400.
     296    d_v(1:klon,1)=-v(1:klon,1)/86400.
     297    ! newtonian relaxation towards temp_newton()
     298    do k=1,klev
     299       temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3
     300       d_t(1:klon,k)=(temp_newton(1:klon,k)-temp(1:klon,k))/5.e5
     301    enddo
     302else
     303   temp_newton(:,:)=0.
     304endif
     305
     306
     307!------------------------------------------------------------
     308! Reevaporation de la pluie
     309!------------------------------------------------------------
     310
     311iflag_ice_thermo=1
     312if ( iflag_reevap == 1 ) then
     313      CALL reevap (klon,klev,iflag_ice_thermo,t_loc,q_loc(:,:,1),q_loc(:,:,2),q_loc(:,:,3), &
     314&                  d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
     315     do k=1,klev
     316        do i=1,klon
     317           t_loc(i,k)=t_loc(i,k)+d_t_eva(i,k)
     318           q_loc(i,k,1)=q_loc(i,k,1)+d_q_eva(i,k)
     319           q_loc(i,k,2)=q_loc(i,k,2)+d_ql_eva(i,k)
     320           q_loc(i,k,3)=q_loc(i,k,3)+d_qi_eva(i,k)
     321           q_loc(i,k,2)=0.
     322           q_loc(i,k,3)=0.
     323        enddo
     324     enddo
     325else
     326     d_t_eva(:,:)=0.
     327     d_q_eva(:,:)=0.
     328     d_ql_eva(:,:)=0.
     329     d_qi_eva(:,:)=0.
     330endif
     331
     332
     333
     334!-----------------------------------------------------------------------
     335! Variables intermédiaires (altitudes, temperature potentielle ...)
     336!-----------------------------------------------------------------------
     337
     338DO k=1,klev
     339   DO i=1,klon
     340      zzlay(i,k)=pphi(i,k)/rg
     341   ENDDO
     342ENDDO
     343DO i=1,klon
     344   zzlev(i,1)=0.
     345ENDDO
     346DO k=2,klev
     347   DO i=1,klon
     348      z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
     349      z2=(paprs(i,k)+pplay(i,k))/(paprs(i,k)-pplay(i,k))
     350      zzlev(i,k)=(z1*zzlay(i,k-1)+z2*zzlay(i,k))/(z1+z2)
     351   ENDDO
     352ENDDO
     353
     354!   Transformation de la temperature en temperature potentielle
     355DO k=1,klev
     356   DO i=1,klon
     357      zpopsk(i,k)=(pplay(i,k)/paprs(i,1))**rkappa
     358      masse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
     359   ENDDO
     360ENDDO
     361DO k=1,klev
     362   DO i=1,klon
     363      h_loc(i,k)=t_loc(i,k)/zpopsk(i,k)
     364      d_h_loc(i,k)=d_t_loc(i,k)/zpopsk(i,k)
     365      d_q_loc(i,k,1)=0.
     366   ENDDO
     367ENDDO
     368
     369!-----------------------------------------------------------------------
     370! Diffusion verticale
     371!-----------------------------------------------------------------------
     372
     373if ( iflag_vdif == 1 ) then
     374      emis(:)=1.
     375      !tsrf=300.
     376      z0m=0.035
     377      z0h=0.035
     378      capcal=1e2
     379      lwrite=.false.
     380      print*,'lwrite ',lwrite
     381            call vdif(klon,klev,                                               &
     382     &                pdtphys,capcal,z0m,z0h,                                  &
     383     &                pplay,paprs,zzlay,zzlev,                                 &
     384     &                u_loc,v_loc,t_loc,h_loc,q_loc,tsrf,emis,                 &
     385     &                d_u_loc,d_v_loc,d_h_loc,d_q_loc,fluxsrf,                 &
     386     &                d_u_vdif,d_v_vdif,d_h_vdif,d_q_vdif,dtsrf,q2,kz_v,kz_h,  &
     387     &                richardson,cdv,cdh,                                      &
     388     &                lwrite)
     389      do k=1,klev
     390         do i=1,klon
     391            d_t_vdif(i,k)=d_h_vdif(i,k)*zpopsk(i,k)
     392            t_loc(i,k)=t_loc(i,k)+d_t_vdif(i,k)*pdtphys
     393            u_loc(i,k)=u_loc(i,k)+d_u_vdif(i,k)*pdtphys
     394            v_loc(i,k)=v_loc(i,k)+d_v_vdif(i,k)*pdtphys
     395            q_loc(i,k,1)=q_loc(i,k,1)+d_q_vdif(i,k,1)*pdtphys
     396         enddo
     397      enddo
     398      do i=1,klon
     399         tsrf(i)=tsrf(i)+dtsrf(i)*pdtphys
     400      enddo
     401else
     402      d_u_vdif(:,:)=0.
     403      d_v_vdif(:,:)=0.
     404      d_t_vdif(:,:)=0.
     405      d_h_vdif(:,:)=0.
     406      d_q_vdif(:,:,1)=0.
     407      kz_v(:,:)=0.
     408      kz_h(:,:)=0.
     409      richardson(:,:)=0.
     410endif
     411
     412!-----------------------------------------------------------------------
     413! Thermiques
     414!-----------------------------------------------------------------------
     415
    119416do k=1,klev
    120   temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3
    121   d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5
     417   do i=1,klon
     418      d_u_the(i,k)=0.
     419      d_v_the(i,k)=0.
     420      d_t_the(i,k)=0.
     421      d_q_the(i,k,1)=0.
     422   enddo
    122423enddo
    123424
     425if ( iflag_thermals > 0 ) then
     426
     427
     428          zqta(:,:)=q_loc(:,:,1)
     429          call caltherm(pdtphys                            &
     430         &      ,pplay,paprs,pphi                          &
     431         &      ,u_loc,v_loc,t_loc,q_loc,debut             &
     432         &      ,f0,zmax0,d_u_the,d_v_the,d_t_the,d_q_the  &
     433         &     ,frac_the,fm_the,zqta,ztv,zpspsk,ztla,zthl  &
     434         &   )
     435   
     436          do k=1,klev
     437             do i=1,klon
     438                t_loc(i,k)=t_loc(i,k)+d_t_the(i,k)
     439                u_loc(i,k)=u_loc(i,k)+d_u_the(i,k)
     440                v_loc(i,k)=v_loc(i,k)+d_v_the(i,k)
     441                q_loc(i,k,1)=q_loc(i,k,1)+d_q_the(i,k,1)
     442             enddo
     443          enddo
     444
     445else
     446          frac_the(:,:)=0.
     447          fm_the(:,:)=0.
     448          ztv(:,:)=t_loc(:,:)
     449          zqta(:,:)=q_loc(:,:,1)
     450          ztla(:,:)=0.
     451          zthl(:,:)=0.
     452          zpspsk(:,:)=(pplay(:,:)/100000.)**RKAPPA
     453 
     454endif
     455
     456!-----------------------------------------------------------------------
     457! Condensation grande échelle
     458!-----------------------------------------------------------------------
     459
     460iflag_cld_th=5
     461ok_ice_sursat=.false.
     462ptconv(:,:)=.false.
     463distcltop=0.
     464temp_cltop=0.
     465beta(:,:)=1.
     466rneb_seri(:,:)=0.
     467do k=1,klev
     468  ratqs(:,k)=ratqsbas+0.5*(ratqshaut-ratqsbas) &
     469  *( tanh( (ratqsp0-pplay(:,k))/ratqsdp) + 1.)
     470enddo
     471
     472
     473if ( iflag_lscp == 1 ) then
     474
     475    call lscp(klon,klev,pdtphys,missing_val,            &
     476     paprs,pplay,t_loc,q_loc,ptconv,ratqs,                      &
     477     d_t_lscp, d_q_lscp, d_ql_lscp, d_qi_lscp, rneb, rneblsvol, rneb_seri,  &
     478     pfraclr,pfracld,                                   &
     479     radocond, radicefrac, rain, snow,                  &
     480     frac_impa, frac_nucl, beta,                        &
     481     prfl, psfl, rhcl, zqta, frac_the,                     &
     482     ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
     483     iflag_ice_thermo, ok_ice_sursat, qsatl, qsats,     &
     484     distcltop,temp_cltop,                               &
     485     qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss,   &
     486     Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
     487     cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
     488
     489
     490          do k=1,klev
     491             do i=1,klon
     492                t_loc(i,k)=t_loc(i,k)+d_t_lscp(i,k)
     493                q_loc(i,k,1)=q_loc(i,k,1)+d_q_lscp(i,k)
     494                q_loc(i,k,2)=q_loc(i,k,2)+d_ql_lscp(i,k)
     495                q_loc(i,k,3)=q_loc(i,k,3)+d_qi_lscp(i,k)
     496             enddo
     497          enddo
     498
     499else
     500     d_t_lscp(:,:)=0.
     501     d_q_lscp(:,:)=0.
     502     d_ql_lscp(:,:)=0.
     503     d_qi_lscp(:,:)=0.
     504     rneb(:,:)=0.
     505     rneblsvol(:,:)=0.
     506     pfraclr(:,:)=0.
     507     pfracld(:,:)=0.
     508     radocond(:,:)=0.
     509     rain(:)=0.
     510     snow(:)=0.
     511     radicefrac(:,:)=0.
     512     rhcl      (:,:)=0.
     513     qsatl     (:,:)=0.
     514     qsats     (:,:)=0.
     515     prfl      (:,:)=0.
     516     psfl      (:,:)=0.
     517     distcltop (:,:)=0.
     518     temp_cltop(:,:)=0.
     519     frac_impa (:,:)=0.
     520     frac_nucl (:,:)=0.
     521     qclr      (:,:)=0.
     522     qcld      (:,:)=0.
     523     qss       (:,:)=0.
     524     qvc       (:,:)=0.
     525     rnebclr   (:,:)=0.
     526     rnebss    (:,:)=0.
     527     gamma_ss  (:,:)=0.
     528     Tcontr    (:,:)=0.
     529     qcontr    (:,:)=0.
     530     qcontr2   (:,:)=0.
     531     fcontrN   (:,:)=0.
     532     fcontrP   (:,:)=0.
     533endif
     534
     535
     536d_u(:,:)=(u_loc(:,:)-u(:,:))/pdtphys
     537d_v(:,:)=(v_loc(:,:)-v(:,:))/pdtphys
     538d_t(:,:)=(t_loc(:,:)-temp(:,:))/pdtphys
     539d_qx(:,:,:)=(q_loc(:,:,:)-qx(:,:,:))/pdtphys
    124540
    125541!------------------------------------------------------------
    … …  
    128544
    129545
    130 call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,t,qx,0.*t,0.*t,0.*t,0.*t,0.*t,0.*t)
     546tsrf_(:)=tsrf(:)
     547if ( iflag_replay == -1 ) then
     548   call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,temp,qx,0.*u,0.*u,0.*u,0.*u,q2,0.*u)
     549else if (iflag_replay == 0 ) then
     550     ! En mode replay, on sort aussi les variables de base
     551! Les lignes qui suivent ont été générées automatiquement avec :
     552! ( for i in `grep -i 'real.*::' physiqex_mod.F90 | sed -e '/^!/d' | grep '(klon,klev' | cut -d: -f3 | cut -d! -f1 | sed -e 's/,/ /g' -e '/rot/d'` ; do echo ' call iophys_ecrit("'$i'",klev,"","",'$i')' ; done ) > physiqex_out.h
     553! ( for i in `grep -i 'real.*::' physiqex_mod.F90 | sed -e '/^!/d' | grep '(klon)' | cut -d: -f3 | cut -d! -f1 | sed -e 's/,/ /g' -e '/rot/d'` ; do echo ' call iophys_ecrit("'$i'",1,"","",'$i')' ; done ) >> physiqex_out.h
     554include "physiqex_out.h"
     555
     556endif
    131557
    132558
    … …  
    136562endif
    137563
     564print*,'Fin physiqex'
    138565
    139566end subroutine physiqex

  • LMDZ6/trunk/libf/phylmd/vdif_kcay.F90

    r4649 r4654  
    22! $Header$
    33
    4 SUBROUTINE vdif_kcay(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, &
     4SUBROUTINE vdif_kcay(klon,klev, dt, g, rconst, plev, temp, zlev, zlay, u, v, &
    55    teta, cd, q2, q2diag, km, kn, ustar, l_mix)
    6   USE dimphy
     6
    77  IMPLICIT NONE
    88
    … …  
    2828
    2929  ! .......................................................................
    30   REAL dt, g, rconst
    31   REAL plev(klon, klev+1), temp(klon, klev)
    32   REAL ustar(klon), snstable
    33   REAL zlev(klon, klev+1)
    34   REAL zlay(klon, klev)
    35   REAL u(klon, klev)
    36   REAL v(klon, klev)
    37   REAL teta(klon, klev)
    38   REAL cd(klon)
    39   REAL q2(klon, klev+1), q2s(klon, klev+1)
    40   REAL q2diag(klon, klev+1)
    41   REAL km(klon, klev+1)
    42   REAL kn(klon, klev+1)
    43   REAL sq(klon), sqz(klon), zz(klon, klev+1), zq, long0(klon)
    44 
    45   INTEGER l_mix, iii
     30  INTEGER, INTENT(IN) :: klon,klev
     31  REAL,INTENT(IN) :: dt, g, rconst
     32  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: plev
     33  REAL,DIMENSION(klon,klev),INTENT(IN) :: temp
     34  REAL,DIMENSION(klon),INTENT(IN) :: ustar(klon)
     35  REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: zlev
     36  REAL,DIMENSION(klon,klev),INTENT(IN) :: zlay
     37  REAL,DIMENSION(klon,klev),INTENT(IN) :: u
     38  REAL,DIMENSION(klon,klev),INTENT(IN) :: v
     39  REAL,DIMENSION(klon,klev),INTENT(IN) :: teta
     40  REAL,DIMENSION(klon),INTENT(IN) :: cd(klon)
     41  REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: q2
     42  REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: q2diag
     43  REAL,DIMENSION(klon,klev+1),INTENT(OUT) :: km
     44  REAL,DIMENSION(klon,klev+1),INTENT(OUT) :: kn
     45  INTEGER, INTENT(OUT) :: l_mix
     46
     47  REAL,DIMENSION(klon) :: sq, sqz, long0
     48  REAL,DIMENSION(klon,klev+1) :: q2s,zz
     49  REAL :: snstable,zq
     50
     51  INTEGER :: iii
    4652  ! .......................................................................
    4753
    … …  
    5662
    5763  ! .......................................................................
    58   INTEGER nlay, nlev, ngrid
    59   REAL unsdz(klon, klev)
    60   REAL unsdzdec(klon, klev+1)
    61   REAL q(klon, klev+1)
     64  INTEGER :: nlay, nlev, ngrid
     65  REAL, DIMENSION(klon,klev) :: unsdz
     66  REAL, DIMENSION(klon, klev+1) :: unsdzdec,q
    6267
    6368  ! .......................................................................
    … …  
    7075
    7176  ! .......................................................................
    72   REAL kmpre(klon, klev+1)
    73   REAL qcstat
    74   REAL q2cstat
    75   REAL sss, sssq
     77  REAL, DIMENSION(klon, klev+1) :: kmpre
     78  REAL :: qcstat
     79  REAL :: q2cstat
     80  REAL :: sss, sssq
    7681  ! .......................................................................
    7782
    … …  
    7984
    8085  ! .......................................................................
    81   REAL long(klon, klev+1)
     86  REAL, DIMENSION(klon, klev+1) :: long
    8287  ! .......................................................................
    8388
    … …  
    97102
    98103  ! .......................................................................
    99   REAL kmq3
    100   REAL kmcstat
    101   REAL knq3
    102   REAL mcstat
    103   REAL m2cstat
    104   REAL m(klon, klev+1)
    105   REAL mpre(klon, klev+1)
    106   REAL m2(klon, klev+1)
    107   REAL n2(klon, klev+1)
     104  REAL :: kmq3
     105  REAL :: kmcstat
     106  REAL :: knq3
     107  REAL :: mcstat
     108  REAL :: m2cstat
     109  REAL, DIMENSION(klon, klev+1) :: m,mpre,m2,n2
    108110  ! .......................................................................
    109111
    … …  
    121123
    122124  ! .......................................................................
    123   REAL gn
    124   REAL gnmin
    125   REAL gnmax
    126   LOGICAL gninf
    127   LOGICAL gnsup
    128   REAL gm
    129   ! REAL ri(klon,klev+1)
    130   REAL sn(klon, klev+1)
    131   REAL snq2(klon, klev+1)
    132   REAL sm(klon, klev+1)
    133   REAL smq2(klon, klev+1)
     125  REAL :: gn,gm
     126  REAL :: gnmin
     127  REAL :: gnmax
     128  LOGICAL :: gninf
     129  LOGICAL :: gnsup
     130  REAL, DIMENSION(klon, klev+1) :: sn, snq2, sm, smq2
    134131  ! .......................................................................
    135132
    … …  
    142139
    143140  ! .......................................................................
    144   REAL kappa
    145   REAL long00
    146   REAL a1, a2, b1, b2, c1
    147   REAL cn1, cn2
    148   REAL cm1, cm2, cm3, cm4
     141  REAL :: kappa
     142  REAL :: long00
     143  REAL :: a1, a2, b1, b2, c1
     144  REAL :: cn1, cn2
     145  REAL :: cm1, cm2, cm3, cm4
    149146  ! .......................................................................
    150147
    … …  
    155152
    156153  ! .......................................................................
    157   REAL termq
    158   REAL termq3
    159   REAL termqm2
    160   REAL termq3m2
     154  REAL :: termq
     155  REAL :: termq3
     156  REAL :: termqm2
     157  REAL :: termq3m2
    161158  ! .......................................................................
    162159
    … …  
    165162
    166163  ! .......................................................................
    167   REAL q2min
    168   REAL q2max
     164  REAL :: q2min
     165  REAL :: q2max
    169166  ! .......................................................................
    170167  ! knmin : borne inferieure de kn
    171168  ! kmmin : borne inferieure de km
    172169  ! .......................................................................
    173   REAL knmin
    174   REAL kmmin
    175   ! .......................................................................
    176   INTEGER ilay, ilev, igrid
    177   REAL tmp1, tmp2
     170  REAL :: knmin
     171  REAL :: kmmin
     172  ! .......................................................................
     173  INTEGER :: ilay, ilev, igrid
     174  REAL :: tmp1, tmp2
    178175  ! .......................................................................
    179176  PARAMETER (kappa=0.4E+0)
    … …  
    202199  PARAMETER (cm4=-9.E+0*a1*a2)
    203200
    204   LOGICAL first
     201  LOGICAL :: first
    205202  SAVE first
    206203  DATA first/.TRUE./
    … …  
    211208
    212209  ! Initialisation de q2
     210ngrid=klon
    213211  nlay = klev
    214212  nlev = klev + 1
    215213
    216   CALL yamada(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, cd, &
    217     q2diag, km, kn, ustar, l_mix)
    218   IF (first .AND. 1==1) THEN
    219     first = .FALSE.
    220     q2 = q2diag
    221   END IF
     214! Initialisation avec un schema d'equilibre
     215! CALL yamada(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, cd, &
     216!   q2diag, km, kn, ustar, l_mix)
     217! IF (first .AND. 1==1) THEN
     218!   first = .FALSE.
     219!   q2 = q2diag
     220! END IF
     221q2diag=0.
    222222
    223223  DO ilev = 2, nlay
Note: See TracChangeset for help on using the changeset viewer.