source: LMDZ6/trunk/libf/phylmd/physiq_mod.F90 @ 4114

Last change on this file since 4114 was 4114, checked in by evignon, 2 years ago

option pour inclure les effets radiatifs des chutes de neige
plus corrections du schema des nuages de phase mixte

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 190.1 KB
RevLine 
[2418]1!
[1279]2! $Id: physiq_mod.F90 4114 2022-03-28 16:23:04Z evignon $
[2418]3!
[1862]4!#define IO_DEBUG
[2418]5MODULE physiq_mod
[766]6
[2469]7  IMPLICIT NONE
[2418]8
9CONTAINS
10
[2469]11  SUBROUTINE physiq (nlon,nlev, &
12       debut,lafin,pdtphys_, &
13       paprs,pplay,pphi,pphis,presnivs, &
14       u,v,rot,t,qx, &
15       flxmass_w, &
16       d_u, d_v, d_t, d_qx, d_ps)
[524]17
[3776]18! For clarity, the "USE" section is now arranged in alphabetical order,
19! with a separate section for CPP keys
20! PLEASE try to follow this rule
21
22    USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
23    USE aero_mod
24    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
25  &      fl_ebil, fl_cor_ebil
[3387]26    USE assert_m, only: assert
[3776]27    USE change_srf_frac_mod
28    USE conf_phys_m, only: conf_phys
29    USE carbon_cycle_mod, ONLY : infocfields_init, RCO2_glo, carbon_cycle_rad
30    USE CFMIP_point_locations   ! IM stations CFMIP
31    USE cmp_seri_mod
32    USE dimphy
33    USE etat0_limit_unstruct_mod
34    USE FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
35    USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
36    USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
[2469]37    USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
38         histwrite, ju2ymds, ymds2ju, getin
[3776]39    USE ioipsl_getin_p_mod, ONLY : getin_p
40    USE indice_sol_mod
[4056]41    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac, nqCO2
[4069]42    USE readTracFiles_mod, ONLY: addPhase
[4056]43    USE strings_mod,  ONLY: strIdx
[3776]44    USE iophy
45    USE limit_read_mod, ONLY : init_limit_read
[3435]46    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid1dTo2d_glo, grid_type, unstructured
[3776]47    USE mod_phys_lmdz_mpi_data, only: is_mpi_root
[2469]48    USE mod_phys_lmdz_para
[3776]49    USE netcdf95, only: nf95_close
50    USE netcdf, only: nf90_fill_real     ! IM for NMC files
51    USE open_climoz_m, only: open_climoz ! ozone climatology from a file
52    USE ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
53    USE pbl_surface_mod, ONLY : pbl_surface
54    USE phyaqua_mod, only: zenang_an
55    USE phystokenc_mod, ONLY: offline, phystokenc
56    USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, &
57         year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour
58!!  USE phys_local_var_mod, ONLY : a long list of variables
59!!              ==> see below, after "CPP Keys" section
60    USE phys_state_var_mod ! Variables sauvegardees de la physique
61    USE phys_output_mod
62    USE phys_output_ctrlout_mod
[3981]63    USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level, &
64         alert_first_call, call_alert, prt_alerte
[3776]65    USE readaerosol_mod, ONLY : init_aero_fromfile
66    USE readaerosolstrato_m, ONLY : init_readaerosolstrato
67    USE radlwsw_m, only: radlwsw
68    USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
69    USE regr_pr_time_av_m, only: regr_pr_time_av
[3792]70    USE surface_data,     ONLY : type_ocean, ok_veget, landice_opt
[3776]71    USE time_phylmdz_mod, only: annee_ref, current_time, day_ini, day_ref, &
72          day_step_phy, itau_phy, pdtphys, raz_date, start_time, update_time
73    USE tracinca_mod, ONLY: config_inca
[2788]74    USE tropopause_m,     ONLY: dyn_tropopause
[4059]75    USE ice_sursat_mod,  ONLY: flight_init, airplane
[3776]76    USE vampir
77    USE VERTICAL_LAYERS_MOD, ONLY: aps,bps, ap, bp
78    USE write_field_phy
[3999]79    USE lscp_mod, ONLY : lscp
[4085]80    USE wake_ini_mod, ONLY : wake_ini
[4089]81    USE thermcell_ini_mod, ONLY : thermcell_ini
[3776]82
83    !USE cmp_seri_mod
84!    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
85!  &      fl_ebil, fl_cor_ebil
86
87!!!!!!!!!!!!!!!!!! "USE" section for CPP keys !!!!!!!!!!!!!!!!!!!!!!!!
88!
89!
[2630]90#ifdef CPP_Dust
[3776]91    USE phytracr_spl_mod, ONLY: phytracr_spl, phytracr_spl_out_init
92    USE phys_output_write_spl_mod
93#else
94    USE phytrac_mod, ONLY : phytrac_init, phytrac
95    USE phys_output_write_mod
[2630]96#endif
[3776]97
98
99#ifdef REPROBUS
100    USE CHEM_REP, ONLY : Init_chem_rep_xjour, &
101         d_q_rep,d_ql_rep,d_qi_rep,ptrop,ttrop, &
102         ztrop, gravit,itroprep, Z1,Z2,fac,B
103#endif
104
105
106#ifdef CPP_RRTM
107    USE YOERAD, ONLY : NRADLP
108    USE YOESW, ONLY : RSUN
109#endif
110
111
[3522]112#ifdef CPP_StratAer
113    USE strataer_mod, ONLY: strataer_init
114#endif
[3776]115
116
117#ifdef CPP_XIOS
118    USE xios, ONLY: xios_update_calendar, xios_context_finalize, &
119            xios_get_field_attr, xios_field_is_active
120    USE wxios, ONLY: missing_val, missing_val_omp
121#endif
122#ifndef CPP_XIOS
123    USE paramLMDZ_phy_mod
124#endif
125!
126!
127!!!!!!!!!!!!!!!!!!  END "USE" for CPP keys !!!!!!!!!!!!!!!!!!!!!!
128
129USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
[2606]130       ! [Variables internes non sauvegardees de la physique]
131       ! Variables locales pour effectuer les appels en serie
[4059]132       t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
[2606]133       ! Dynamic tendencies (diagnostics)
[4059]134       d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn,d_rneb_dyn, &
[2606]135       d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, &
136       ! Physic tendencies
137       d_t_con,d_q_con,d_u_con,d_v_con, &
138       d_tr, &                              !! to be removed?? (jyg)
139       d_t_wake,d_q_wake, &
140       d_t_lwr,d_t_lw0,d_t_swr,d_t_sw0, &
141       d_t_ajsb,d_q_ajsb, &
142       d_t_ajs,d_q_ajs,d_u_ajs,d_v_ajs, &
143       d_t_ajs_w,d_q_ajs_w, &
144       d_t_ajs_x,d_q_ajs_x, &
145       !
[2705]146       d_t_eva,d_q_eva,d_ql_eva,d_qi_eva, &
[2606]147       d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc, &
148       d_t_lscst,d_q_lscst, &
149       d_t_lscth,d_q_lscth, &
150       plul_st,plul_th, &
151       !
152       d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
[4056]153       d_t_vdf_x, d_t_vdf_w, &
154       d_q_vdf_x, d_q_vdf_w, &
[2606]155       d_ts, &
156       !
157       d_t_oli,d_u_oli,d_v_oli, &
158       d_t_oro,d_u_oro,d_v_oro, &
[2897]159       d_t_oro_gw,d_u_oro_gw,d_v_oro_gw, &
[2606]160       d_t_lif,d_u_lif,d_v_lif, &
161       d_t_ec, &
162       !
163       du_gwd_hines,dv_gwd_hines,d_t_hin, &
164       dv_gwd_rando,dv_gwd_front, &
165       east_gwstress,west_gwstress, &
166       d_q_ch4, &
167       !  Special RRTM
168       ZLWFT0_i,ZSWFT0_i,ZFLDN0,  &
169       ZFLUP0,ZFSDN0,ZFSUP0,      &
170       !
171       topswad_aero,solswad_aero,   &
172       topswai_aero,solswai_aero,   &
173       topswad0_aero,solswad0_aero, &
174       !LW additional
175       toplwad_aero,sollwad_aero,   &
176       toplwai_aero,sollwai_aero,   &
177       toplwad0_aero,sollwad0_aero, &
178       !
179       topsw_aero,solsw_aero,       &
180       topsw0_aero,solsw0_aero,     &
181       topswcf_aero,solswcf_aero,   &
182       tausum_aero,tau3d_aero,      &
[2854]183       drytausum_aero,              &
[2606]184       !
185       !variables CFMIP2/CMIP5
186       topswad_aerop, solswad_aerop,   &
187       topswai_aerop, solswai_aerop,   &
188       topswad0_aerop, solswad0_aerop, &
189       topsw_aerop, topsw0_aerop,      &
190       solsw_aerop, solsw0_aerop,      &
191       topswcf_aerop, solswcf_aerop,   &
192       !LW diagnostics
193       toplwad_aerop, sollwad_aerop,   &
194       toplwai_aerop, sollwai_aerop,   &
195       toplwad0_aerop, sollwad0_aerop, &
196       !
197       ptstar, pt0, slp, &
198       !
199       bils, &
200       !
201       cldh, cldl,cldm, cldq, cldt,      &
202       JrNt,                             &
203       dthmin, evap, fder, plcl, plfc,   &
204       prw, prlw, prsw,                  &
205       s_lcl, s_pblh, s_pblt, s_therm,   &
206       cdragm, cdragh,                   &
207       zustar, zu10m, zv10m, rh2m, qsat2m, &
[3817]208       zq2m, zt2m, zn2mout, weak_inversion, &
[2606]209       zt2m_min_mon, zt2m_max_mon,   &         ! pour calcul_divers.h
210       t2m_min_mon, t2m_max_mon,  &            ! pour calcul_divers.h
211       !
212       s_pblh_x, s_pblh_w, &
213       s_lcl_x, s_lcl_w,   &
214       !
215       slab_wfbils, tpot, tpote,               &
216       ue, uq, ve, vq, zxffonte,               &
[3257]217       uwat, vwat,                             &
[2606]218       zxfqcalving, zxfluxlat,                 &
219       zxrunofflic,                            &
220       zxtsol, snow_lsc, zxfqfonte, zxqsurf,   &
[3888]221       delta_qsurf,                            &
[2606]222       rain_lsc, rain_num,                     &
223       !
224       sens_x, sens_w, &
225       zxfluxlat_x, zxfluxlat_w, &
226       !
[4009]227       pbl_tke_input, tke_dissip, l_mix, wprime,&
[2606]228       t_therm, q_therm, u_therm, v_therm, &
229       cdragh_x, cdragh_w, &
230       cdragm_x, cdragm_w, &
231       kh, kh_x, kh_w, &
232       !
[2730]233       wake_k, &
[3080]234       alp_wake, &
[2635]235       wake_h, wake_omg, &
236                       ! tendencies of delta T and delta q:
237       d_deltat_wk, d_deltaq_wk, &         ! due to wakes
238       d_deltat_wk_gw, d_deltaq_wk_gw, &   ! due to wake induced gravity waves
239       d_deltat_vdf, d_deltaq_vdf, &       ! due to vertical diffusion
240       d_deltat_the, d_deltaq_the, &       ! due to thermals
241       d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection
242                       ! tendencies of wake fractional area and wake number per unit area:
[3208]243       d_s_wk,  d_dens_a_wk,  d_dens_wk, &  ! due to wakes
244!!!       d_s_vdf, d_dens_a_vdf, d_dens_vdf, & ! due to vertical diffusion
245!!!       d_s_the, d_dens_a_the, d_dens_the, & ! due to thermals
[2635]246       !                                 
[3148]247       ptconv, ratqsc, &
[2824]248       wbeff, convoccur, zmax_th, &
[2606]249       sens, flwp, fiwp,  &
[3080]250       alp_bl_conv,alp_bl_det,  &
[2606]251       alp_bl_fluct_m,alp_bl_fluct_tke,  &
252       alp_bl_stat, n2, s2,  &
253       proba_notrig, random_notrig,  &
[3956]254!!       cv_gen,  &  !moved to phys_state_var_mod
[2606]255       !
[3134]256       dnwd0,  &
257       omega,  &
[2606]258       epmax_diag,  &
[3134]259       !    Deep convective variables used in phytrac
260       pmflxr, pmflxs,  &
[3496]261       wdtrainA, wdtrainS, wdtrainM,  &
[3134]262       upwd, dnwd, &
[2606]263       ep,  &
[3134]264       da, mp, &
265       phi, &
266       wght_cvfd, &
267       phi2, &
268       d1a, dam, &
269       ev, &
270       elij, &
[3496]271       qtaa, &
[3134]272       clw, &
273       epmlmMm, eplaMm, &
274       sij, &
[3387]275       !
[2606]276       cldemi,  &
277       cldfra, cldtau, fiwc,  &
278       fl, re, flwc,  &
279       ref_liq, ref_ice, theta,  &
280       ref_liq_pi, ref_ice_pi,  &
[3780]281       zphi, zx_rh, zx_rhl, zx_rhi,  &
[2606]282       pmfd, pmfu,  &
283       !
284       t2m, fluxlat,  &
285       fsollw, evap_pot,  &
286       fsolsw, wfbils, wfbilo,  &
[2670]287       wfevap, wfrain, wfsnow,  & 
[3134]288       prfl, psfl, fraca, Vprecip,  &
[2606]289       zw2,  &
[3387]290       !
[2606]291       fluxu, fluxv,  &
292       fluxt,  &
[3387]293       !
[2606]294       uwriteSTD, vwriteSTD, &                !pour calcul_STDlev.h
295       wwriteSTD, phiwriteSTD, &              !pour calcul_STDlev.h
296       qwriteSTD, twriteSTD, rhwriteSTD, &    !pour calcul_STDlev.h
[3387]297       !
[2606]298       beta_prec,  &
299       rneb,  &
[2968]300       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
[2606]301       !
[2630]302
[782]303
[3387]304    IMPLICIT NONE
[2469]305    !>======================================================================
306    !!
307    !! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
308    !!
309    !! Objet: Moniteur general de la physique du modele
310    !!AA      Modifications quant aux traceurs :
311    !!AA                  -  uniformisation des parametrisations ds phytrac
312    !!AA                  -  stockage des moyennes des champs necessaires
313    !!AA                     en mode traceur off-line
314    !!======================================================================
315    !!   CLEFS CPP POUR LES IO
316    !!   =====================
[1352]317#define histNMC
[2469]318    !!======================================================================
319    !!    modif   ( P. Le Van ,  12/10/98 )
320    !!
321    !!  Arguments:
322    !!
323    !! nlon----input-I-nombre de points horizontaux
324    !! nlev----input-I-nombre de couches verticales, doit etre egale a klev
325    !! debut---input-L-variable logique indiquant le premier passage
326    !! lafin---input-L-variable logique indiquant le dernier passage
327    !! jD_cur       -R-jour courant a l'appel de la physique (jour julien)
328    !! jH_cur       -R-heure courante a l'appel de la physique (jour julien)
329    !! pdtphys-input-R-pas d'integration pour la physique (seconde)
330    !! paprs---input-R-pression pour chaque inter-couche (en Pa)
331    !! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
332    !! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
333    !! pphis---input-R-geopotentiel du sol
334    !! presnivs-input_R_pressions approximat. des milieux couches ( en PA)
335    !! u-------input-R-vitesse dans la direction X (de O a E) en m/s
336    !! v-------input-R-vitesse Y (de S a N) en m/s
337    !! t-------input-R-temperature (K)
338    !! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
339    !! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
340    !! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
[2496]341    !! d_ql_dyn-input-R-tendance dynamique pour "ql" (kg/kg/s)
342    !! d_qs_dyn-input-R-tendance dynamique pour "qs" (kg/kg/s)
[2469]343    !! flxmass_w -input-R- flux de masse verticale
344    !! d_u-----output-R-tendance physique de "u" (m/s/s)
345    !! d_v-----output-R-tendance physique de "v" (m/s/s)
346    !! d_t-----output-R-tendance physique de "t" (K/s)
347    !! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
348    !! d_ps----output-R-tendance physique de la pression au sol
349    !!======================================================================
350    integer jjmp1
351    !  parameter (jjmp1=jjm+1-1/jjm) ! => (jjmp1=nbp_lat-1/(nbp_lat-1))
352    !  integer iip1
353    !  parameter (iip1=iim+1)
[782]354
[2469]355    include "regdim.h"
356    include "dimsoil.h"
357    include "clesphys.h"
[4089]358    include "alpale.h"
[3011]359    include "dimpft.h"
[2469]360    !======================================================================
[3479]361    LOGICAL, SAVE :: ok_volcan ! pour activer les diagnostics volcaniques
[3600]362    !$OMP THREADPRIVATE(ok_volcan)
[4056]363    INTEGER, SAVE :: flag_volc_surfstrat ! pour imposer le cool/heat rate à la surf/strato
[3989]364    !$OMP THREADPRIVATE(flag_volc_surfstrat)
[2469]365    LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE
366    PARAMETER (ok_cvl=.TRUE.)
367    LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
368    PARAMETER (ok_gust=.FALSE.)
[3600]369    INTEGER, SAVE :: iflag_radia     ! active ou non le rayonnement (MPL)
[2469]370    !$OMP THREADPRIVATE(iflag_radia)
371    !======================================================================
372    LOGICAL check ! Verifier la conservation du modele en eau
373    PARAMETER (check=.FALSE.)
374    LOGICAL ok_stratus ! Ajouter artificiellement les stratus
375    PARAMETER (ok_stratus=.FALSE.)
376    !======================================================================
377    REAL amn, amx
378    INTEGER igout
379    !======================================================================
[3317]380    ! Clef iflag_cycle_diurne controlant l'activation du cycle diurne:
[2469]381    ! en attente du codage des cles par Fred
[3317]382    ! iflag_cycle_diurne est initialise par conf_phys et se trouve
383    ! dans clesphys.h (IM)
[2469]384    !======================================================================
385    ! Modele thermique du sol, a activer pour le cycle diurne:
386    !cc      LOGICAL soil_model
387    !cc      PARAMETER (soil_model=.FALSE.)
388    !======================================================================
389    ! Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
390    ! le calcul du rayonnement est celle apres la precipitation des nuages.
391    ! Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
392    ! la condensation et la precipitation. Cette cle augmente les impacts
393    ! radiatifs des nuages.
394    !cc      LOGICAL new_oliq
395    !cc      PARAMETER (new_oliq=.FALSE.)
396    !======================================================================
397    ! Clefs controlant deux parametrisations de l'orographie:
398    !c      LOGICAL ok_orodr
399    !cc      PARAMETER (ok_orodr=.FALSE.)
400    !cc      LOGICAL ok_orolf
401    !cc      PARAMETER (ok_orolf=.FALSE.)
402    !======================================================================
403    LOGICAL ok_journe ! sortir le fichier journalier
[3461]404    SAVE ok_journe
[2469]405    !$OMP THREADPRIVATE(ok_journe)
406    !
407    LOGICAL ok_mensuel ! sortir le fichier mensuel
[3461]408    SAVE ok_mensuel
[2469]409    !$OMP THREADPRIVATE(ok_mensuel)
410    !
411    LOGICAL ok_instan ! sortir le fichier instantane
[3461]412    SAVE ok_instan
[2469]413    !$OMP THREADPRIVATE(ok_instan)
414    !
415    LOGICAL ok_LES ! sortir le fichier LES
[3461]416    SAVE ok_LES                           
[2469]417    !$OMP THREADPRIVATE(ok_LES)                 
418    !
419    LOGICAL callstats ! sortir le fichier stats
[3461]420    SAVE callstats                           
[2469]421    !$OMP THREADPRIVATE(callstats)                 
422    !
423    LOGICAL ok_region ! sortir le fichier regional
424    PARAMETER (ok_region=.FALSE.)
425    !======================================================================
[3461]426    REAL seuil_inversion
427    SAVE seuil_inversion
[2469]428    !$OMP THREADPRIVATE(seuil_inversion)
[3461]429    INTEGER iflag_ratqs
430    SAVE iflag_ratqs
[2469]431    !$OMP THREADPRIVATE(iflag_ratqs)
432    real facteur
[1507]433
[2469]434    REAL wmax_th(klon)
435    REAL tau_overturning_th(klon)
[878]436
[3461]437    INTEGER lmax_th(klon)
438    INTEGER limbas(klon)
439    REAL ratqscth(klon,klev)
440    REAL ratqsdiff(klon,klev)
441    REAL zqsatth(klon,klev)
[878]442
[2469]443    !======================================================================
444    !
445    INTEGER ivap          ! indice de traceurs pour vapeur d'eau
446    PARAMETER (ivap=1)
447    INTEGER iliq          ! indice de traceurs pour eau liquide
448    PARAMETER (iliq=2)
449    INTEGER isol          ! indice de traceurs pour eau glace
450    PARAMETER (isol=3)
[4059]451    INTEGER irneb         ! indice de traceurs pour fraction nuageuse LS (optional)
452    PARAMETER (irneb=4)   
[2469]453    !
454    !
455    ! Variables argument:
456    !
457    INTEGER nlon
458    INTEGER nlev
459    REAL,INTENT(IN) :: pdtphys_
460    ! NB: pdtphys to be used in physics is in time_phylmdz_mod
461    LOGICAL debut, lafin
462    REAL paprs(klon,klev+1)
463    REAL pplay(klon,klev)
464    REAL pphi(klon,klev)
465    REAL pphis(klon)
466    REAL presnivs(klev)
[2799]467!JLD    REAL znivsig(klev)
468!JLD    real pir
[719]469
[2469]470    REAL u(klon,klev)
471    REAL v(klon,klev)
[2333]472
[2469]473    REAL, intent(in):: rot(klon, klev)
474    ! relative vorticity, in s-1, needed for frontal waves
[2333]475
[2469]476    REAL t(klon,klev),thetal(klon,klev)
477    ! thetal: ligne suivante a decommenter si vous avez les fichiers
478    !     MPL 20130625
479    ! fth_fonctions.F90 et parkind1.F90
480    ! sinon thetal=theta
481    !     REAL fth_thetae,fth_thetav,fth_thetal
482    REAL qx(klon,klev,nqtot)
483    REAL flxmass_w(klon,klev)
484    REAL d_u(klon,klev)
485    REAL d_v(klon,klev)
486    REAL d_t(klon,klev)
487    REAL d_qx(klon,klev,nqtot)
488    REAL d_ps(klon)
[2897]489  ! variables pour tend_to_tke
490    REAL duadd(klon,klev)
491    REAL dvadd(klon,klev)
492    REAL dtadd(klon,klev)
493
[2271]494#ifndef CPP_XIOS
[2997]495    REAL, SAVE :: missing_val=nf90_fill_real
[2271]496#endif
[3134]497!!   Variables moved to phys_local_var_mod
498!!    ! Variables pour le transport convectif
499!!    real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
500!!    real wght_cvfd(klon,klev)
501!!    ! Variables pour le lessivage convectif
502!!    ! RomP >>>
503!!    real phi2(klon,klev,klev)
504!!    real d1a(klon,klev),dam(klon,klev)
505!!    real ev(klon,klev)
506!!    real clw(klon,klev),elij(klon,klev,klev)
507!!    real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
508!!    ! RomP <<<
[2469]509    !IM definition dynamique o_trac dans phys_output_open
510    !      type(ctrl_out) :: o_trac(nqtot)
[524]511
[2469]512    ! variables a une pression donnee
513    !
514    include "declare_STDlev.h"
515    !
516    !
517    include "radopt.h"
518    !
519    !
520    INTEGER debug
521    INTEGER n
522    !ym      INTEGER npoints
523    !ym      PARAMETER(npoints=klon)
524    !
525    INTEGER nregISCtot
526    PARAMETER(nregISCtot=1)
527    !
528    ! imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties
529    ! sur 1 region rectangulaire y compris pour 1 point
530    ! imin_debut : indice minimum de i; nbpti : nombre de points en
531    ! direction i (longitude)
532    ! jmin_debut : indice minimum de j; nbptj : nombre de points en
533    ! direction j (latitude)
[2799]534!JLD    INTEGER imin_debut, nbpti
535!JLD    INTEGER jmin_debut, nbptj
[2469]536    !IM: region='3d' <==> sorties en global
537    CHARACTER*3 region
538    PARAMETER(region='3d')
[3461]539    LOGICAL ok_hf
[2469]540    !
[3461]541    SAVE ok_hf
[2469]542    !$OMP THREADPRIVATE(ok_hf)
[524]543
[3461]544    INTEGER, PARAMETER :: longcles=20
545    REAL, SAVE :: clesphy0(longcles)
[2469]546    !$OMP THREADPRIVATE(clesphy0)
547    !
548    ! Variables propres a la physique
[3461]549    INTEGER, SAVE :: itap         ! compteur pour la physique
[2469]550    !$OMP THREADPRIVATE(itap)
[2235]551
[2469]552    INTEGER, SAVE :: abortphy=0   ! Reprere si on doit arreter en fin de phys
553    !$OMP THREADPRIVATE(abortphy)
554    !
[3461]555    REAL,SAVE ::  solarlong0
[2469]556    !$OMP THREADPRIVATE(solarlong0)
[987]557
[2469]558    !
559    !  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
560    !
561    !IM 141004     REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
562    REAL zulow(klon),zvlow(klon)
563    !
564    INTEGER igwd,idx(klon),itest(klon)
565    !
566    !      REAL,allocatable,save :: run_off_lic_0(:)
567    ! !$OMP THREADPRIVATE(run_off_lic_0)
568    !ym      SAVE run_off_lic_0
569    !KE43
570    ! Variables liees a la convection de K. Emanuel (sb):
571    !
[3461]572    REAL, SAVE :: bas, top             ! cloud base and top levels
[2469]573    !$OMP THREADPRIVATE(bas, top)
574    !------------------------------------------------------------------
575    ! Upmost level reached by deep convection and related variable (jyg)
576    !
577    INTEGER izero
578    INTEGER k_upper_cv
579    !------------------------------------------------------------------
[3153]580    ! Compteur de l'occurence de cvpas=1
581    INTEGER Ncvpaseq1
582    SAVE Ncvpaseq1
583    !$OMP THREADPRIVATE(Ncvpaseq1)
[2469]584    !
585    !==========================================================================
586    !CR04.12.07: on ajoute les nouvelles variables du nouveau schema
587    !de convection avec poches froides
588    ! Variables li\'ees \`a la poche froide (jyg)
[879]589
[3496]590!!    REAL mipsh(klon,klev)  ! mass flux shed by the adiab ascent at each level
591!!      Moved to phys_state_var_mod
[2469]592    !
593    REAL wape_prescr, fip_prescr
594    INTEGER it_wape_prescr
595    SAVE wape_prescr, fip_prescr, it_wape_prescr
596    !$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
597    !
598    ! variables supplementaires de concvl
599    REAL Tconv(klon,klev)
[3134]600!!    variable moved to phys_local_var_mod
601!!    REAL sij(klon,klev,klev)
[2812]602!!    !
603!!    ! variables pour tester la conservation de l'energie dans concvl
604!!    REAL, DIMENSION(klon,klev)     :: d_t_con_sat
605!!    REAL, DIMENSION(klon,klev)     :: d_q_con_sat
606!!    REAL, DIMENSION(klon,klev)     :: dql_sat
[970]607
[3461]608    REAL, SAVE :: alp_bl_prescr=0.
609    REAL, SAVE :: ale_bl_prescr=0.
610    REAL, SAVE :: wake_s_min_lsp=0.1
[2469]611    !$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
612    !$OMP THREADPRIVATE(wake_s_min_lsp)
[970]613
[3461]614    REAL ok_wk_lsp(klon)
[1516]615
[2469]616    !RC
617    ! Variables li\'ees \`a la poche froide (jyg et rr)
[879]618
[2635]619    INTEGER,  SAVE               :: iflag_wake_tend  ! wake: if =0, then wake state variables are
620                                                     ! updated within calwake
621    !$OMP THREADPRIVATE(iflag_wake_tend)
[3000]622    INTEGER,  SAVE               :: iflag_alp_wk_cond=0 ! wake: if =0, then Alp_wk is the average lifting
623                                                        ! power provided by the wakes; else, Alp_wk is the
624                                                        ! lifting power conditionned on the presence of a
625                                                        ! gust-front in the grid cell.
626    !$OMP THREADPRIVATE(iflag_alp_wk_cond)
[3489]627
[2635]628    REAL t_w(klon,klev),q_w(klon,klev) ! temperature and moisture profiles in the wake region
629    REAL t_x(klon,klev),q_x(klon,klev) ! temperature and moisture profiles in the off-wake region
[879]630
[2469]631    REAL wake_dth(klon,klev)        ! wake : temp pot difference
[879]632
[2469]633    REAL wake_omgbdth(klon,klev)    ! Wake : flux of Delta_Theta
634    ! transported by LS omega
635    REAL wake_dp_omgb(klon,klev)    ! Wake : vertical gradient of
636    ! large scale omega
637    REAL wake_dtKE(klon,klev)       ! Wake : differential heating
638    ! (wake - unpertubed) CONV
639    REAL wake_dqKE(klon,klev)       ! Wake : differential moistening
640    ! (wake - unpertubed) CONV
641    REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
642    REAL wake_spread(klon,klev)     ! spreading term in wake_delt
643    !
644    !pourquoi y'a pas de save??
645    !
[2730]646!!!    INTEGER, SAVE, DIMENSION(klon)   :: wake_k
647!!!    !$OMP THREADPRIVATE(wake_k)
[2469]648    !
649    !jyg<
650    !cc      REAL wake_pe(klon)              ! Wake potential energy - WAPE
651    !>jyg
[879]652
[3000]653    REAL wake_fip_0(klon)           ! Average Front Incoming Power (unconditionned)
[2469]654    REAL wake_gfl(klon)             ! Gust Front Length
[2635]655!!!    REAL wake_dens(klon)         ! moved to phys_state_var_mod
[2469]656    !
657    !
658    REAL dt_dwn(klon,klev)
659    REAL dq_dwn(klon,klev)
660    REAL M_dwn(klon,klev)
661    REAL M_up(klon,klev)
662    REAL dt_a(klon,klev)
663    REAL dq_a(klon,klev)
664    REAL d_t_adjwk(klon,klev)                !jyg
665    REAL d_q_adjwk(klon,klev)                !jyg
666    LOGICAL,SAVE :: ok_adjwk=.FALSE.
667    !$OMP THREADPRIVATE(ok_adjwk)
[2882]668    INTEGER,SAVE :: iflag_adjwk=0            !jyg
669    !$OMP THREADPRIVATE(iflag_adjwk)         !jyg
[2657]670    REAL,SAVE :: oliqmax=999.,oicemax=999.
671    !$OMP THREADPRIVATE(oliqmax,oicemax)
[2469]672    REAL, SAVE :: alp_offset
673    !$OMP THREADPRIVATE(alp_offset)
[3150]674    REAL, SAVE :: dtcon_multistep_max=1.e6
675    !$OMP THREADPRIVATE(dtcon_multistep_max)
676    REAL, SAVE :: dqcon_multistep_max=1.e6
677    !$OMP THREADPRIVATE(dqcon_multistep_max)
678
[2897]679 
[2469]680    !
681    !RR:fin declarations poches froides
682    !==========================================================================
[1032]683
[2469]684    REAL ztv(klon,klev),ztva(klon,klev)
685    REAL zpspsk(klon,klev)
686    REAL ztla(klon,klev),zqla(klon,klev)
687    REAL zthl(klon,klev)
[1638]688
[2469]689    !cc nrlmd le 10/04/2012
[1638]690
[2469]691    !--------Stochastic Boundary Layer Triggering: ALE_BL--------
692    !---Propri\'et\'es du thermiques au LCL
693    real zlcl_th(klon)          ! Altitude du LCL calcul\'e
694    ! continument (pcon dans
695    ! thermcell_main.F90)
696    real fraca0(klon)           ! Fraction des thermiques au LCL
697    real w0(klon)               ! Vitesse des thermiques au LCL
698    real w_conv(klon)           ! Vitesse verticale de grande \'echelle au LCL
699    real tke0(klon,klev+1)      ! TKE au d\'ebut du pas de temps
700    real therm_tke_max0(klon)   ! TKE dans les thermiques au LCL
701    real env_tke_max0(klon)     ! TKE dans l'environnement au LCL
[1638]702
[2799]703!JLD    !---D\'eclenchement stochastique
704!JLD    integer :: tau_trig(klon)
[1638]705
[2469]706    REAL,SAVE :: random_notrig_max=1.
707    !$OMP THREADPRIVATE(random_notrig_max)
[2294]708
[2469]709    !--------Statistical Boundary Layer Closure: ALP_BL--------
710    !---Profils de TKE dans et hors du thermique
711    real therm_tke_max(klon,klev)   ! Profil de TKE dans les thermiques
712    real env_tke_max(klon,klev)     ! Profil de TKE dans l'environnement
[1638]713
[2897]714    !-------Activer les tendances de TKE due a l'orograp??ie---------
715     INTEGER, SAVE :: addtkeoro
716    !$OMP THREADPRIVATE(addtkeoro)
717     REAL, SAVE :: alphatkeoro
718    !$OMP THREADPRIVATE(alphatkeoro)
719     LOGICAL, SAVE :: smallscales_tkeoro
720    !$OMP THREADPRIVATE(smallscales_tkeoro)
[1638]721
[2897]722
723
[2469]724    !cc fin nrlmd le 10/04/2012
[782]725
[2469]726    ! Variables locales pour la couche limite (al1):
727    !
728    !Al1      REAL pblh(klon)           ! Hauteur de couche limite
729    !Al1      SAVE pblh
730    !34EK
731    !
732    ! Variables locales:
733    !
734    !AA
735    !AA  Pour phytrac
736    REAL u1(klon)             ! vents dans la premiere couche U
737    REAL v1(klon)             ! vents dans la premiere couche V
[524]738
[2469]739    !@$$      LOGICAL offline           ! Controle du stockage ds "physique"
740    !@$$      PARAMETER (offline=.false.)
741    !@$$      INTEGER physid
742    REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
743    REAL frac_nucl(klon,klev) ! idem (nucleation)
744    ! RomP >>>
745    REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
746    ! RomP <<<
[2832]747    REAL          :: calday
[2068]748
[2469]749    !IM cf FH pour Tiedtke 080604
750    REAL rain_tiedtke(klon),snow_tiedtke(klon)
751    !
752    !IM 050204 END
753    REAL devap(klon) ! evaporation et sa derivee
754    REAL dsens(klon) ! chaleur sensible et sa derivee
[1279]755
[2469]756    !
757    ! Conditions aux limites
758    !
759    !
760    REAL :: day_since_equinox
761    ! Date de l'equinoxe de printemps
762    INTEGER, parameter :: mth_eq=3, day_eq=21
763    REAL :: jD_eq
[1279]764
[3461]765    LOGICAL, parameter :: new_orbit = .TRUE.
[524]766
[2469]767    !
768    INTEGER lmt_pas
769    SAVE lmt_pas                ! frequence de mise a jour
770    !$OMP THREADPRIVATE(lmt_pas)
771    real zmasse(klon, nbp_lev),exner(klon, nbp_lev)
772    !     (column-density of mass of air in a cell, in kg m-2)
773    real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
[1797]774
[2469]775    !IM sorties
776    REAL un_jour
777    PARAMETER(un_jour=86400.)
778    INTEGER itapm1 !pas de temps de la physique du(es) mois precedents
779    SAVE itapm1    !mis a jour le dernier pas de temps du mois en cours
780    !$OMP THREADPRIVATE(itapm1)
781    !======================================================================
782    !
783    ! Declaration des procedures appelees
784    !
785    EXTERNAL angle     ! calculer angle zenithal du soleil
786    EXTERNAL alboc     ! calculer l'albedo sur ocean
787    EXTERNAL ajsec     ! ajustement sec
788    EXTERNAL conlmd    ! convection (schema LMD)
789    !KE43
790    EXTERNAL conema3  ! convect4.3
791    EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
792    !AA
793    ! JBM (3/14) fisrtilp_tr not loaded
794    ! EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
795    !                          ! stockage des coefficients necessaires au
796    !                          ! lessivage OFF-LINE et ON-LINE
797    EXTERNAL hgardfou  ! verifier les temperatures
798    EXTERNAL nuage     ! calculer les proprietes radiatives
799    !C      EXTERNAL o3cm      ! initialiser l'ozone
800    EXTERNAL orbite    ! calculer l'orbite terrestre
801    EXTERNAL phyetat0  ! lire l'etat initial de la physique
802    EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
803    EXTERNAL suphel    ! initialiser certaines constantes
804    EXTERNAL transp    ! transport total de l'eau et de l'energie
805    !IM
806    EXTERNAL haut2bas  !variables de haut en bas
807    EXTERNAL ini_undefSTD  !initialise a 0 une variable a 1 niveau de pression
808    EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
809    !     EXTERNAL moy_undefSTD  !moyenne d'1 var a 1 niveau de pression
810    ! EXTERNAL moyglo_aire
811    ! moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
812    ! par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
813    !
814    !
815    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
816    ! Local variables
817    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
818    !
819    REAL rhcl(klon,klev)    ! humiditi relative ciel clair
820    REAL dialiq(klon,klev)  ! eau liquide nuageuse
821    REAL diafra(klon,klev)  ! fraction nuageuse
822    REAL cldliq(klon,klev)  ! eau liquide nuageuse
823    !
824    !XXX PB
825    REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
826    !
827    REAL zxfluxt(klon, klev)
828    REAL zxfluxq(klon, klev)
829    REAL zxfluxu(klon, klev)
830    REAL zxfluxv(klon, klev)
[1797]831
[2469]832    ! Le rayonnement n'est pas calcule tous les pas, il faut donc
833    !                      sauvegarder les sorties du rayonnement
834    !ym      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
835    !ym      SAVE  sollwdownclr, toplwdown, toplwdownclr
836    !ym      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
837    !
838    INTEGER itaprad
839    SAVE itaprad
840    !$OMP THREADPRIVATE(itaprad)
841    !
842    REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
843    REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
844    !
[2799]845#ifdef INCA
[2469]846    REAL zxsnow_dummy(klon)
[2799]847#endif
[2469]848    REAL zsav_tsol(klon)
849    !
850    REAL dist, rmu0(klon), fract(klon)
851    REAL zrmu0(klon), zfract(klon)
852    REAL zdtime, zdtime1, zdtime2, zlongi
853    !
854    REAL qcheck
855    REAL z_avant(klon), z_apres(klon), z_factor(klon)
856    LOGICAL zx_ajustq
857    !
[2799]858    REAL za
859    REAL zx_t, zx_qs, zdelta, zcor
[2469]860    real zqsat(klon,klev)
861    !
[4056]862    INTEGER i, k, iq, j, nsrf, ll, l, itr
[2469]863    !
864    REAL t_coup
865    PARAMETER (t_coup=234.0)
[1797]866
[2469]867    !ym A voir plus tard !!
868    !ym      REAL zx_relief(iim,jjmp1)
869    !ym      REAL zx_aire(iim,jjmp1)
870    !
871    ! Grandeurs de sorties
872    REAL s_capCL(klon)
873    REAL s_oliqCL(klon), s_cteiCL(klon)
874    REAL s_trmb1(klon), s_trmb2(klon)
875    REAL s_trmb3(klon)
[2707]876
877    ! La convection n'est pas calculee tous les pas, il faut donc
878    !                      sauvegarder les sorties de la convection
879    !ym      SAVE 
880    !ym      SAVE 
881    !ym      SAVE 
882    !
[2730]883    INTEGER itapcv, itapwk
884    SAVE itapcv, itapwk
885    !$OMP THREADPRIVATE(itapcv, itapwk)
[2707]886
[2469]887    !KE43
888    ! Variables locales pour la convection de K. Emanuel (sb):
[524]889
[2469]890    REAL tvp(klon,klev)       ! virtual temp of lifted parcel
891    CHARACTER*40 capemaxcels  !max(CAPE)
[1412]892
[2469]893    REAL rflag(klon)          ! flag fonctionnement de convect
894    INTEGER iflagctrl(klon)          ! flag fonctionnement de convect
[904]895
[2469]896    ! -- convect43:
897    INTEGER ntra              ! nb traceurs pour convect4.3
898    REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
899    REAL dplcldt(klon), dplcldr(klon)
900    !?     .     condm_con(klon,klev),conda_con(klon,klev),
901    !?     .     mr_con(klon,klev),ep_con(klon,klev)
902    !?     .    ,sadiab(klon,klev),wadiab(klon,klev)
903    ! --
904    !34EK
905    !
906    ! Variables du changement
907    !
908    ! con: convection
909    ! lsc: condensation a grande echelle (Large-Scale-Condensation)
910    ! ajs: ajustement sec
911    ! eva: evaporation de l'eau liquide nuageuse
912    ! vdf: couche limite (Vertical DiFfusion)
[2611]913    !
[2469]914    ! tendance nulles
[2812]915    REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
916    REAL, dimension(klon)     :: dsig0, ddens0
917    INTEGER, dimension(klon)  :: wkoccur1
[2801]918    ! tendance buffer pour appel de add_phys_tend
919    REAL, DIMENSION(klon,klev)  :: d_q_ch4_dtime
[2611]920    !
921    ! Flag pour pouvoir ne pas ajouter les tendances.
922    ! Par defaut, les tendances doivente etre ajoutees et
923    ! flag_inhib_tend = 0
924    ! flag_inhib_tend > 0 : tendances non ajoutees, avec un nombre
925    ! croissant de print quand la valeur du flag augmente
926    !!! attention, ce flag doit etre change avec prudence !!!
927    INTEGER :: flag_inhib_tend = 0 !  0 is the default value
928!!    INTEGER :: flag_inhib_tend = 2
[3134]929    !
930    ! Logical switch to a bug : reseting to 0 convective variables at the
931    ! begining of physiq.
932    LOGICAL, SAVE :: ok_bug_cv_trac = .TRUE.
933    !$OMP THREADPRIVATE(ok_bug_cv_trac)
[3180]934    !
935    ! Logical switch to a bug : changing wake_deltat when thermals are active
936    ! even when there are no wakes.
937    LOGICAL, SAVE :: ok_bug_split_th = .TRUE.
938    !$OMP THREADPRIVATE(ok_bug_split_th)
[524]939
[2469]940    !
941    !********************************************************
942    !     declarations
[524]943
[2469]944    !********************************************************
945    !IM 081204 END
946    !
947    REAL pen_u(klon,klev), pen_d(klon,klev)
948    REAL pde_u(klon,klev), pde_d(klon,klev)
949    INTEGER kcbot(klon), kctop(klon), kdtop(klon)
950    !
[3461]951    REAL ratqsbas,ratqshaut,tau_ratqs
952    SAVE ratqsbas,ratqshaut,tau_ratqs
[2469]953    !$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
[2534]954    REAL, SAVE :: ratqsp0=50000., ratqsdp=20000.
955    !$OMP THREADPRIVATE(ratqsp0, ratqsdp)
[644]956
[2469]957    ! Parametres lies au nouveau schema de nuages (SB, PDF)
[3461]958    REAL, SAVE :: fact_cldcon
959    REAL, SAVE :: facttemps
960    !$OMP THREADPRIVATE(fact_cldcon,facttemps)
961    LOGICAL, SAVE :: ok_newmicro
[2469]962    !$OMP THREADPRIVATE(ok_newmicro)
[524]963
[3461]964    INTEGER, SAVE :: iflag_cld_th
[2469]965    !$OMP THREADPRIVATE(iflag_cld_th)
[2877]966!IM logical ptconv(klon,klev)  !passe dans phys_local_var_mod
[2469]967    !IM cf. AM 081204 BEG
[3461]968    LOGICAL ptconvth(klon,klev)
[3999]969
970    REAL picefra(klon,klev)
[2469]971    !IM cf. AM 081204 END
972    !
973    ! Variables liees a l'ecriture de la bande histoire physique
974    !
975    !======================================================================
976    !
977    !
[2799]978!JLD    integer itau_w   ! pas de temps ecriture = itap + itau_phy
[2469]979    !
980    !
981    ! Variables locales pour effectuer les appels en serie
982    !
983    !IM RH a 2m (la surface)
984    REAL Lheat
[524]985
[2469]986    INTEGER        length
987    PARAMETER    ( length = 100 )
988    REAL tabcntr0( length       )
989    !
[2799]990!JLD    INTEGER ndex2d(nbp_lon*nbp_lat)
[2469]991    !IM
992    !
993    !IM AMIP2 BEG
[2799]994!JLD    REAL moyglo, mountor
[2469]995    !IM 141004 BEG
996    REAL zustrdr(klon), zvstrdr(klon)
997    REAL zustrli(klon), zvstrli(klon)
998    REAL zustrph(klon), zvstrph(klon)
999    REAL aam, torsfc
1000    !IM 141004 END
1001    !IM 190504 BEG
1002    !  INTEGER imp1jmp1
1003    !  PARAMETER(imp1jmp1=(iim+1)*jjmp1)
1004    !ym A voir plus tard
1005    !  REAL zx_tmp((nbp_lon+1)*nbp_lat)
1006    !  REAL airedyn(nbp_lon+1,nbp_lat)
1007    !IM 190504 END
[2799]1008!JLD    LOGICAL ok_msk
1009!JLD    REAL msk(klon)
[2469]1010    !ym A voir plus tard
1011    !ym      REAL zm_wo(jjmp1, klev)
1012    !IM AMIP2 END
1013    !
1014    REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
1015    REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
[2799]1016!JLD    REAL zx_tmp_2d(nbp_lon,nbp_lat)
1017!JLD    REAL zx_lon(nbp_lon,nbp_lat)
1018!JLD    REAL zx_lat(nbp_lon,nbp_lat)
[2469]1019    !
[2630]1020    INTEGER nid_ctesGCM
1021    SAVE nid_ctesGCM
1022    !$OMP THREADPRIVATE(nid_ctesGCM)
[2469]1023    !
1024    !IM 280405 BEG
1025    !  INTEGER nid_bilKPins, nid_bilKPave
1026    !  SAVE nid_bilKPins, nid_bilKPave
1027    !  !$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
1028    !
1029    REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
1030    REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
1031    REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
1032    REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
1033    !
[2799]1034!JLD    REAL zjulian
1035!JLD    SAVE zjulian
1036!JLD!$OMP THREADPRIVATE(zjulian)
[2590]1037
[2799]1038!JLD    INTEGER nhori, nvert
1039!JLD    REAL zsto
1040!JLD    REAL zstophy, zout
[2068]1041
[3981]1042    CHARACTER (LEN=20) :: modname='physiq_mod'
[4056]1043    CHARACTER*80 abort_message
[3461]1044    LOGICAL, SAVE ::  ok_sync, ok_sync_omp
[2469]1045    !$OMP THREADPRIVATE(ok_sync)
[3461]1046    REAL date0
[524]1047
[2469]1048    ! essai writephys
[3461]1049    INTEGER fid_day, fid_mth, fid_ins
1050    PARAMETER (fid_ins = 1, fid_day = 2, fid_mth = 3)
1051    INTEGER prof2d_on, prof3d_on, prof2d_av, prof3d_av
1052    PARAMETER (prof2d_on = 1, prof3d_on = 2, prof2d_av = 3, prof3d_av = 4)
[2469]1053    REAL ztsol(klon)
1054    REAL q2m(klon,nbsrf)  ! humidite a 2m
[524]1055
[2469]1056    !IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
1057    CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
[2799]1058    CHARACTER*40 tinst, tave
[2469]1059    REAL cldtaupi(klon,klev) ! Cloud optical thickness for
1060    ! pre-industrial (pi) aerosols
[524]1061
[2863]1062    INTEGER :: naero
[2469]1063    ! Aerosol optical properties
1064    CHARACTER*4, DIMENSION(naero_grp) :: rfname
1065    REAL, DIMENSION(klon,klev)     :: mass_solu_aero ! total mass
1066    ! concentration
1067    ! for all soluble
1068    ! aerosols[ug/m3]
1069    REAL, DIMENSION(klon,klev)     :: mass_solu_aero_pi
1070    ! - " - (pre-industrial value)
[1279]1071
[2469]1072    ! Parameters
1073    LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not
[2738]1074    LOGICAL ok_alw            ! Apply aerosol LW effect or not
[2469]1075    LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013)
1076    REAL bl95_b0, bl95_b1   ! Parameter in Boucher and Lohmann (1995)
[2738]1077    SAVE ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1
1078    !$OMP THREADPRIVATE(ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1)
[2469]1079    LOGICAL, SAVE :: aerosol_couple ! true  : calcul des aerosols dans INCA
1080    ! false : lecture des aerosol dans un fichier
1081    !$OMP THREADPRIVATE(aerosol_couple)   
[3338]1082    LOGICAL, SAVE :: chemistry_couple ! true  : use INCA chemistry O3
1083    ! false : use offline chemistry O3
1084    !$OMP THREADPRIVATE(chemistry_couple)   
[2469]1085    INTEGER, SAVE :: flag_aerosol
1086    !$OMP THREADPRIVATE(flag_aerosol)
[2644]1087    LOGICAL, SAVE :: flag_bc_internal_mixture
1088    !$OMP THREADPRIVATE(flag_bc_internal_mixture)
[2469]1089    !
1090    !--STRAT AEROSOL
[2530]1091    INTEGER, SAVE :: flag_aerosol_strat
[2469]1092    !$OMP THREADPRIVATE(flag_aerosol_strat)
[3412]1093    !
1094    !--INTERACTIVE AEROSOL FEEDBACK ON RADIATION
1095    LOGICAL, SAVE :: flag_aer_feedback
1096    !$OMP THREADPRIVATE(flag_aer_feedback)
1097
[2469]1098    !c-fin STRAT AEROSOL
1099    !
1100    ! Declaration des constantes et des fonctions thermodynamiques
1101    !
[3461]1102    LOGICAL,SAVE :: first=.TRUE.
[2469]1103    !$OMP THREADPRIVATE(first)
[1279]1104
[2788]1105    ! VARIABLES RELATED TO OZONE CLIMATOLOGIES ; all are OpenMP shared
1106    ! Note that pressure vectors are in Pa and in stricly ascending order
1107    INTEGER,SAVE :: read_climoz                ! Read ozone climatology
[2469]1108    !     (let it keep the default OpenMP shared attribute)
1109    !     Allowed values are 0, 1 and 2
1110    !     0: do not read an ozone climatology
1111    !     1: read a single ozone climatology that will be used day and night
1112    !     2: read two ozone climatologies, the average day and night
1113    !     climatology and the daylight climatology
[2788]1114    INTEGER,SAVE :: ncid_climoz                ! NetCDF file identifier
1115    REAL, POINTER, SAVE :: press_cen_climoz(:) ! Pressure levels
1116    REAL, POINTER, SAVE :: press_edg_climoz(:) ! Edges of pressure intervals
1117    REAL, POINTER, SAVE :: time_climoz(:)      ! Time vector
1118    CHARACTER(LEN=13), PARAMETER :: vars_climoz(2) &
1119                                  = ["tro3         ","tro3_daylight"]
1120    ! vars_climoz(1:read_climoz): variables names in climoz file.
[2819]1121    ! vars_climoz(1:read_climoz-2) if read_climoz>2 (temporary)
1122    REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for
1123                 ! the ozone fields, old method.
[1279]1124
[2469]1125    include "YOMCST.h"
1126    include "YOETHF.h"
1127    include "FCTTRE.h"
1128    !IM 100106 BEG : pouvoir sortir les ctes de la physique
1129    include "conema3.h"
1130    include "fisrtilp.h"
1131    include "nuage.h"
1132    include "compbl.h"
1133    !IM 100106 END : pouvoir sortir les ctes de la physique
1134    !
1135    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1136    ! Declarations pour Simulateur COSP
1137    !============================================================
[3511]1138    real :: mr_ozone(klon,klev), phicosp(klon,klev)
[3370]1139
[2469]1140    !IM stations CFMIP
1141    INTEGER, SAVE :: nCFMIP
1142    !$OMP THREADPRIVATE(nCFMIP)
1143    INTEGER, PARAMETER :: npCFMIP=120
1144    INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
1145    REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
1146    !$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
1147    INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
1148    REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
1149    !$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
1150    INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
1151    !$OMP THREADPRIVATE(iGCM, jGCM)
1152    logical, dimension(nfiles)            :: phys_out_filestations
1153    logical, parameter :: lNMC=.FALSE.
[1539]1154
[2469]1155    !IM betaCRF
1156    REAL, SAVE :: pfree, beta_pbl, beta_free
1157    !$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
1158    REAL, SAVE :: lon1_beta,  lon2_beta, lat1_beta, lat2_beta
1159    !$OMP THREADPRIVATE(lon1_beta,  lon2_beta, lat1_beta, lat2_beta)
1160    LOGICAL, SAVE :: mskocean_beta
1161    !$OMP THREADPRIVATE(mskocean_beta)
1162    REAL, dimension(klon, klev) :: beta ! facteur sur cldtaurad et
1163    ! cldemirad pour evaluer les
1164    ! retros liees aux CRF
1165    REAL, dimension(klon, klev) :: cldtaurad   ! epaisseur optique
1166    ! pour radlwsw pour
1167    ! tester "CRF off"
1168    REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique
1169    ! pour radlwsw pour
1170    ! tester "CRF off"
1171    REAL, dimension(klon, klev) :: cldemirad   ! emissivite pour
1172    ! radlwsw pour tester
1173    ! "CRF off"
1174    REAL, dimension(klon, klev) :: cldfrarad   ! fraction nuageuse
[1735]1175
[3418]1176#ifdef INCA
1177    ! set de variables utilisees pour l'initialisation des valeurs provenant de INCA
1178    REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_tauinca
1179    REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_pizinca
1180    REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_cginca
1181    REAL, DIMENSION(klon,klev,nbands) :: init_ccminca
1182#endif
1183    REAL, DIMENSION(klon,nbtr) :: init_source
1184
[3048]1185    !lwoff=y : offset LW CRE for radiation code and other schemes
1186    REAL, SAVE :: betalwoff
1187    !OMP THREADPRIVATE(betalwoff)
1188!
[2469]1189    INTEGER :: nbtr_tmp ! Number of tracer inside concvl
1190    REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
[2784]1191    REAL, dimension(klon,klev) :: ch_in ! Condensed humidity entering in phytrac (eau liquide)
[2469]1192    integer iostat
[1539]1193
[4009]1194    REAL, dimension(klon,klev+1) :: tke_dissip_ave, l_mix_ave, wprime_ave
[2469]1195    REAL zzz
1196    !albedo SB >>>
[3461]1197    REAL,DIMENSION(6), SAVE :: SFRWL
[3435]1198!$OMP THREADPRIVATE(SFRWL)
[2469]1199    !albedo SB <<<
[1955]1200
[2485]1201    !--OB variables for mass fixer (hard coded for now)
[3461]1202    LOGICAL, PARAMETER :: mass_fixer=.FALSE.
1203    REAL qql1(klon),qql2(klon),corrqql
[2476]1204
[3110]1205    REAL pi
1206
1207    pi = 4. * ATAN(1.)
1208
[3981]1209    ! set-up call to alerte function
1210    call_alert = (alert_first_call .AND. is_master)
1211   
[2469]1212    ! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
1213    jjmp1=nbp_lat
[2344]1214
[2469]1215    !======================================================================
1216    ! Gestion calendrier : mise a jour du module phys_cal_mod
1217    !
1218    pdtphys=pdtphys_
1219    CALL update_time(pdtphys)
[3435]1220    phys_tstep=NINT(pdtphys)
1221#ifdef CPP_XIOS
1222    IF (.NOT. debut .AND. is_omp_master) CALL xios_update_calendar(itap+1)
1223#endif
[1355]1224
[2469]1225    !======================================================================
1226    ! Ecriture eventuelle d'un profil verticale en entree de la physique.
1227    ! Utilise notamment en 1D mais peut etre active egalement en 3D
1228    ! en imposant la valeur de igout.
1229    !======================================================================d
[2692]1230    IF (prt_level.ge.1) THEN
[2469]1231       igout=klon/2+1/klon
1232       write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1233       write(lunout,*) 'igout, lat, lon ',igout, latitude_deg(igout), &
1234            longitude_deg(igout)
1235       write(lunout,*) &
1236            'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1237       write(lunout,*) &
1238            nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
[879]1239
[2469]1240       write(lunout,*) 'paprs, play, phi, u, v, t'
[2692]1241       DO k=1,klev
[2469]1242          write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
1243               u(igout,k),v(igout,k),t(igout,k)
[2692]1244       ENDDO
[2469]1245       write(lunout,*) 'ovap (g/kg),  oliq (g/kg)'
[2692]1246       DO k=1,klev
[2469]1247          write(lunout,*) qx(igout,k,1)*1000,qx(igout,k,2)*1000.
[2692]1248       ENDDO
1249    ENDIF
[879]1250
[2769]1251    ! Quick check on pressure levels:
[3461]1252    CALL assert(paprs(:, nbp_lev + 1) < paprs(:, nbp_lev), &
[2769]1253            "physiq_mod paprs bad order")
[879]1254
[2692]1255    IF (first) THEN
[3435]1256       CALL init_etat0_limit_unstruct
1257       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
[2469]1258       !CR:nvelles variables convection/poches froides
1259
[3461]1260       WRITE(lunout,*) '================================================='
1261       WRITE(lunout,*) 'Allocation des variables locales et sauvegardees'
1262       WRITE(lunout,*) '================================================='
[2692]1263       CALL phys_local_var_init
[2469]1264       !
1265       !     appel a la lecture du run.def physique
[2692]1266       CALL conf_phys(ok_journe, ok_mensuel, &
[2469]1267            ok_instan, ok_hf, &
1268            ok_LES, &
1269            callstats, &
1270            solarlong0,seuil_inversion, &
1271            fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
1272            iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
[3989]1273            ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat, aerosol_couple, &
1274            chemistry_couple, flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
[2644]1275            flag_bc_internal_mixture, bl95_b0, bl95_b1, &
[2469]1276                                ! nv flags pour la convection et les
1277                                ! poches froides
1278            read_climoz, &
1279            alp_offset)
[2692]1280       CALL phys_state_var_init(read_climoz)
1281       CALL phys_output_var_init
[3522]1282       IF (read_climoz>=1 .AND. create_etat0_limit .AND. grid_type==unstructured) &
1283          CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
[3435]1284
[3522]1285#ifdef CPP_StratAer
1286       CALL strataer_init
1287#endif
1288
[4059]1289       !!CALL flight_init
1290
[2469]1291       print*, '================================================='
1292       !
1293       !CR: check sur le nb de traceurs de l eau
[2692]1294       IF ((iflag_ice_thermo.gt.0).and.(nqo==2)) THEN
[2469]1295          WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
1296               '(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
[3531]1297          abort_message='see above'
1298          CALL abort_physic(modname,abort_message,1)
[2692]1299       ENDIF
[2224]1300
[4062]1301       IF (ok_ice_sursat.AND.(iflag_ice_thermo.EQ.0)) THEN
1302          WRITE (lunout, *) ' ok_ice_sursat=y requires iflag_ice_thermo=1 as well'
[4059]1303          abort_message='see above'
1304          CALL abort_physic(modname,abort_message,1)
1305       ENDIF
1306
[4062]1307       IF (ok_ice_sursat.AND.(nqo.NE.4)) THEN
1308          WRITE (lunout, *) ' ok_ice_sursat=y requires 4 H2O tracers ', &
[4098]1309               '(H2Ov, H2Ol, H2Oi, H2Or) but nqo=', nqo, '. Might as well stop here.'
[4059]1310          abort_message='see above'
1311          CALL abort_physic(modname,abort_message,1)
1312       ENDIF
1313
[4062]1314       IF (ok_plane_h2o.AND..NOT.ok_ice_sursat) THEN
1315          WRITE (lunout, *) ' ok_plane_h2o=y requires ok_ice_sursat=y '
1316          abort_message='see above'
1317          CALL abort_physic(modname,abort_message,1)
1318       ENDIF
1319
1320       IF (ok_plane_contrail.AND..NOT.ok_ice_sursat) THEN
1321          WRITE (lunout, *) ' ok_plane_contrail=y requires ok_ice_sursat=y '
1322          abort_message='see above'
1323          CALL abort_physic(modname,abort_message,1)
1324       ENDIF
1325
[3154]1326       Ncvpaseq1 = 0
[2469]1327       dnwd0=0.0
1328       ftd=0.0
1329       fqd=0.0
1330       cin=0.
1331       !ym Attention pbase pas initialise dans concvl !!!!
1332       pbase=0
1333       !IM 180608
[904]1334
[2469]1335       itau_con=0
[3461]1336       first=.FALSE.
[1797]1337
[2692]1338    ENDIF  ! first
[1797]1339
[2469]1340    !ym => necessaire pour iflag_con != 2   
1341    pmfd(:,:) = 0.
1342    pen_u(:,:) = 0.
1343    pen_d(:,:) = 0.
1344    pde_d(:,:) = 0.
1345    pde_u(:,:) = 0.
1346    aam=0.
1347    d_t_adjwk(:,:)=0
1348    d_q_adjwk(:,:)=0
[1797]1349
[2469]1350    alp_bl_conv(:)=0.
[2245]1351
[2469]1352    torsfc=0.
1353    forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
[1797]1354
[644]1355
[2469]1356    IF (debut) THEN
1357       CALL suphel ! initialiser constantes et parametres phys.
[3327]1358! tau_gl : constante de rappel de la temperature a la surface de la glace - en
1359       tau_gl=5.
1360       CALL getin_p('tau_gl', tau_gl)
1361! tau_gl : constante de rappel de la temperature a la surface de la glace - en
1362! secondes
1363       tau_gl=86400.*tau_gl
[3461]1364       WRITE(lunout,*) 'debut physiq_mod tau_gl=',tau_gl
1365
[3000]1366       CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond)
[2469]1367       CALL getin_p('random_notrig_max',random_notrig_max)
[2882]1368       CALL getin_p('ok_adjwk',ok_adjwk)
1369       IF (ok_adjwk) iflag_adjwk=2  ! for compatibility with older versions
1370       ! iflag_adjwk: ! 0 = Default: no convective adjustment of w-region
1371                      ! 1 => convective adjustment but state variables are unchanged
1372                      ! 2 => convective adjustment and state variables are changed
1373       CALL getin_p('iflag_adjwk',iflag_adjwk)
[3150]1374       CALL getin_p('dtcon_multistep_max',dtcon_multistep_max)
1375       CALL getin_p('dqcon_multistep_max',dqcon_multistep_max)
[2613]1376       CALL getin_p('oliqmax',oliqmax)
[2657]1377       CALL getin_p('oicemax',oicemax)
[2534]1378       CALL getin_p('ratqsp0',ratqsp0)
1379       CALL getin_p('ratqsdp',ratqsdp)
[2635]1380       iflag_wake_tend = 0
1381       CALL getin_p('iflag_wake_tend',iflag_wake_tend)
[2799]1382       ok_bad_ecmwf_thermo=.TRUE. ! By default thermodynamical constants are set
1383                                  ! in rrtm/suphec.F90 (and rvtmp2 is set to 0).
1384       CALL getin_p('ok_bad_ecmwf_thermo',ok_bad_ecmwf_thermo)
[3134]1385       CALL getin_p('ok_bug_cv_trac',ok_bug_cv_trac)
[3180]1386       CALL getin_p('ok_bug_split_th',ok_bug_split_th)
[2799]1387       fl_ebil = 0 ! by default, conservation diagnostics are desactivated
1388       CALL getin_p('fl_ebil',fl_ebil)
1389       fl_cor_ebil = 0 ! by default, no correction to ensure energy conservation
1390       CALL getin_p('fl_cor_ebil',fl_cor_ebil)
[2984]1391       iflag_phytrac = 1 ! by default we do want to call phytrac
[2973]1392       CALL getin_p('iflag_phytrac',iflag_phytrac)
[3776]1393#ifdef CPP_Dust
1394       IF (iflag_phytrac.EQ.0) THEN
1395         WRITE(lunout,*) 'In order to run with SPLA, iflag_phytrac will be forced to 1'
1396         iflag_phytrac = 1
1397       ENDIF
[4056]1398#endif
[3011]1399       nvm_lmdz = 13
1400       CALL getin_p('NVM',nvm_lmdz)
[3387]1401
[3461]1402       WRITE(lunout,*) 'iflag_alp_wk_cond=',  iflag_alp_wk_cond
1403       WRITE(lunout,*) 'random_ntrig_max=',   random_notrig_max
1404       WRITE(lunout,*) 'ok_adjwk=',           ok_adjwk
1405       WRITE(lunout,*) 'iflag_adjwk=',        iflag_adjwk
1406       WRITE(lunout,*) 'qtcon_multistep_max=',dtcon_multistep_max
1407       WRITE(lunout,*) 'qdcon_multistep_max=',dqcon_multistep_max
1408       WRITE(lunout,*) 'ratqsp0=',            ratqsp0
1409       WRITE(lunout,*) 'ratqsdp=',            ratqsdp
1410       WRITE(lunout,*) 'iflag_wake_tend=',    iflag_wake_tend
1411       WRITE(lunout,*) 'ok_bad_ecmwf_thermo=',ok_bad_ecmwf_thermo
1412       WRITE(lunout,*) 'ok_bug_cv_trac=',     ok_bug_cv_trac
1413       WRITE(lunout,*) 'ok_bug_split_th=',    ok_bug_split_th
1414       WRITE(lunout,*) 'fl_ebil=',            fl_ebil
1415       WRITE(lunout,*) 'fl_cor_ebil=',        fl_cor_ebil
1416       WRITE(lunout,*) 'iflag_phytrac=',      iflag_phytrac
1417       WRITE(lunout,*) 'NVM=',                nvm_lmdz
1418
[3387]1419       !--PC: defining fields to be exchanged between LMDz, ORCHIDEE and NEMO
1420       WRITE(lunout,*) 'Call to infocfields from physiq'
1421       CALL infocfields_init
1422
[2469]1423    ENDIF
[878]1424
[2692]1425    IF (prt_level.ge.1) print *,'CONVERGENCE PHYSIQUE THERM 1 '
[1279]1426
[2469]1427    !======================================================================
1428    ! Gestion calendrier : mise a jour du module phys_cal_mod
1429    !
1430    !     CALL phys_cal_update(jD_cur,jH_cur)
[1279]1431
[2469]1432    !
1433    ! Si c'est le debut, il faut initialiser plusieurs choses
1434    !          ********
1435    !
1436    IF (debut) THEN
1437       !rv CRinitialisation de wght_th et lalim_conv pour la
1438       !definition de la couche alimentation de la convection a partir
1439       !des caracteristiques du thermique
1440       wght_th(:,:)=1.
1441       lalim_conv(:)=1
1442       !RC
1443       ustar(:,:)=0.
[2569]1444!       u10m(:,:)=0.
1445!       v10m(:,:)=0.
[2469]1446       rain_con(:)=0.
1447       snow_con(:)=0.
1448       topswai(:)=0.
1449       topswad(:)=0.
1450       solswai(:)=0.
1451       solswad(:)=0.
[959]1452
[2469]1453       wmax_th(:)=0.
1454       tau_overturning_th(:)=0.
[645]1455
[3865]1456       IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN
[2469]1457          ! jg : initialisation jusqu'au ces variables sont dans restart
1458          ccm(:,:,:) = 0.
1459          tau_aero(:,:,:,:) = 0.
1460          piz_aero(:,:,:,:) = 0.
1461          cg_aero(:,:,:,:) = 0.
[2372]1462
[2469]1463          config_inca='none' ! default
1464          CALL getin_p('config_inca',config_inca)
[2372]1465
[2469]1466       ELSE
1467          config_inca='none' ! default
[2692]1468       ENDIF
[782]1469
[3435]1470       tau_aero(:,:,:,:) = 1.e-15
1471       piz_aero(:,:,:,:) = 1.
1472       cg_aero(:,:,:,:)  = 0.
1473
[2469]1474       IF (aerosol_couple .AND. (config_inca /= "aero" &
1475            .AND. config_inca /= "aeNP ")) THEN
1476          abort_message &
1477               = 'if aerosol_couple is activated, config_inca need to be ' &
1478               // 'aero or aeNP'
1479          CALL abort_physic (modname,abort_message,1)
1480       ENDIF
[2443]1481
[2469]1482       rnebcon0(:,:) = 0.0
1483       clwcon0(:,:) = 0.0
1484       rnebcon(:,:) = 0.0
1485       clwcon(:,:) = 0.0
[1863]1486
[2469]1487       !
1488       print*,'iflag_coupl,iflag_clos,iflag_wake', &
1489            iflag_coupl,iflag_clos,iflag_wake
[3317]1490       print*,'iflag_cycle_diurne', iflag_cycle_diurne
[2469]1491       !
1492       IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
1493          abort_message = 'Tiedtke needs iflag_cld_th=-2 or -1'
1494          CALL abort_physic (modname,abort_message,1)
1495       ENDIF
1496       !
1497       !
1498       ! Initialiser les compteurs:
1499       !
1500       itap    = 0
1501       itaprad = 0
[2707]1502       itapcv = 0
[2730]1503       itapwk = 0
[878]1504
[2469]1505       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1506       !! Un petit travail \`a faire ici.
1507       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[878]1508
[2692]1509       IF (iflag_pbl>1) THEN
[2469]1510          PRINT*, "Using method MELLOR&YAMADA"
[2692]1511       ENDIF
[956]1512
[2469]1513       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1514       ! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans
1515       ! phylmd plutot que dyn3d
1516       ! Attention : la version precedente n'etait pas tres propre.
1517       ! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1518       ! pour obtenir le meme resultat.
[2731]1519!jyg for fh<
[3435]1520       WRITE(lunout,*) 'Pas de temps phys_tstep pdtphys ',phys_tstep,pdtphys
1521       IF (abs(phys_tstep-pdtphys)>1.e-10) THEN
[2731]1522          abort_message='pas de temps doit etre entier en seconde pour orchidee et XIOS'
1523          CALL abort_physic(modname,abort_message,1)
1524       ENDIF
1525!>jyg
[3435]1526       IF (MOD(NINT(86400./phys_tstep),nbapp_rad).EQ.0) THEN
1527          radpas = NINT( 86400./phys_tstep)/nbapp_rad
[2469]1528       ELSE
1529          WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1530               'multiple de nbapp_rad'
1531          WRITE(lunout,*) 'changer nbapp_rad ou alors commenter ce test ', &
1532               'mais 1+1<>2'
1533          abort_message='nbre de pas de temps physique n est pas multiple ' &
1534               // 'de nbapp_rad'
[2692]1535          CALL abort_physic(modname,abort_message,1)
[2469]1536       ENDIF
[3435]1537       IF (nbapp_cv .EQ. 0) nbapp_cv=86400./phys_tstep
1538       IF (nbapp_wk .EQ. 0) nbapp_wk=86400./phys_tstep
[2730]1539       print *,'physiq, nbapp_cv, nbapp_wk ',nbapp_cv,nbapp_wk
[3435]1540       IF (MOD(NINT(86400./phys_tstep),nbapp_cv).EQ.0) THEN
1541          cvpas_0 = NINT( 86400./phys_tstep)/nbapp_cv
[3150]1542          cvpas = cvpas_0
[2707]1543       print *,'physiq, cvpas ',cvpas
1544       ELSE
1545          WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1546               'multiple de nbapp_cv'
1547          WRITE(lunout,*) 'changer nbapp_cv ou alors commenter ce test ', &
1548               'mais 1+1<>2'
1549          abort_message='nbre de pas de temps physique n est pas multiple ' &
1550               // 'de nbapp_cv'
[3461]1551          CALL abort_physic(modname,abort_message,1)
[2707]1552       ENDIF
[3435]1553       IF (MOD(NINT(86400./phys_tstep),nbapp_wk).EQ.0) THEN
1554          wkpas = NINT( 86400./phys_tstep)/nbapp_wk
[3457]1555!       print *,'physiq, wkpas ',wkpas
[2730]1556       ELSE
1557          WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1558               'multiple de nbapp_wk'
1559          WRITE(lunout,*) 'changer nbapp_wk ou alors commenter ce test ', &
1560               'mais 1+1<>2'
1561          abort_message='nbre de pas de temps physique n est pas multiple ' &
1562               // 'de nbapp_wk'
[3461]1563          CALL abort_physic(modname,abort_message,1)
[2730]1564       ENDIF
[2469]1565       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3465]1566       CALL init_iophy_new(latitude_deg,longitude_deg)
[524]1567
[3435]1568          !===================================================================
1569          !IM stations CFMIP
1570          nCFMIP=npCFMIP
1571          OPEN(98,file='npCFMIP_param.data',status='old', &
1572               form='formatted',iostat=iostat)
1573          IF (iostat == 0) THEN
1574             READ(98,*,end=998) nCFMIP
1575998          CONTINUE
1576             CLOSE(98)
1577             CONTINUE
1578             IF(nCFMIP.GT.npCFMIP) THEN
1579                print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1580                CALL abort_physic("physiq", "", 1)
1581             ELSE
1582                print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1583             ENDIF
1584
1585             !
1586             ALLOCATE(tabCFMIP(nCFMIP))
1587             ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1588             ALLOCATE(tabijGCM(nCFMIP))
1589             ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1590             ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1591             !
1592             ! lecture des nCFMIP stations CFMIP, de leur numero
1593             ! et des coordonnees geographiques lonCFMIP, latCFMIP
1594             !
1595             CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP,  &
1596                  lonCFMIP, latCFMIP)
1597             !
1598             ! identification des
1599             ! 1) coordonnees lonGCM, latGCM des points CFMIP dans la
1600             ! grille de LMDZ
1601             ! 2) indices points tabijGCM de la grille physique 1d sur
1602             ! klon points
1603             ! 3) indices iGCM, jGCM de la grille physique 2d
1604             !
1605             CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1606                  tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1607             !
1608          ELSE
1609             ALLOCATE(tabijGCM(0))
1610             ALLOCATE(lonGCM(0), latGCM(0))
1611             ALLOCATE(iGCM(0), jGCM(0))
1612          ENDIF
1613
1614#ifdef CPP_IOIPSL
1615
1616       !$OMP MASTER
1617       ! FH : if ok_sync=.true. , the time axis is written at each time step
1618       ! in the output files. Only at the end in the opposite case
[3461]1619       ok_sync_omp=.FALSE.
[3435]1620       CALL getin('ok_sync',ok_sync_omp)
1621       CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
1622            iGCM,jGCM,lonGCM,latGCM, &
1623            jjmp1,nlevSTD,clevSTD,rlevSTD, phys_tstep,ok_veget, &
1624            type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
1625            ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1626            read_climoz, phys_out_filestations, &
[3630]1627            aerosol_couple, &
[3435]1628            flag_aerosol_strat, pdtphys, paprs, pphis,  &
1629            pplay, lmax_th, ptconv, ptconvth, ivap,  &
1630            d_u, d_t, qx, d_qx, zmasse, ok_sync_omp)
1631       !$OMP END MASTER
1632       !$OMP BARRIER
1633       ok_sync=ok_sync_omp
1634
1635       freq_outNMC(1) = ecrit_files(7)
1636       freq_outNMC(2) = ecrit_files(8)
1637       freq_outNMC(3) = ecrit_files(9)
1638       WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1639       WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1640       WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1641
1642#ifndef CPP_XIOS
1643       CALL ini_paramLMDZ_phy(phys_tstep,nid_ctesGCM)
1644#endif
1645
1646#endif
1647       ecrit_reg = ecrit_reg * un_jour
1648       ecrit_tra = ecrit_tra * un_jour
1649
1650       !XXXPB Positionner date0 pour initialisation de ORCHIDEE
1651       date0 = jD_ref
1652       WRITE(*,*) 'physiq date0 : ',date0
1653       !
1654
1655!       CALL create_climoz(read_climoz)
[3472]1656      IF (.NOT. create_etat0_limit) CALL init_aero_fromfile(flag_aerosol)  !! initialise aero from file for XIOS interpolation (unstructured_grid)
1657      IF (.NOT. create_etat0_limit) CALL init_readaerosolstrato(flag_aerosol_strat)  !! initialise aero strato from file for XIOS interpolation (unstructured_grid)
[3435]1658
1659#ifdef CPP_COSP
1660      IF (ok_cosp) THEN
[3726]1661!           DO k = 1, klev
1662!             DO i = 1, klon
1663!               phicosp(i,k) = pphi(i,k) + pphis(i)
1664!             ENDDO
1665!           ENDDO
[3465]1666        CALL phys_cosp(itap,phys_tstep,freq_cosp, &
[3435]1667               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1668               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1669               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1670               JrNt,ref_liq,ref_ice, &
1671               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1672               zu10m,zv10m,pphis, &
1673               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1674               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1675               prfl(:,1:klev),psfl(:,1:klev), &
1676               pmflxr(:,1:klev),pmflxs(:,1:klev), &
1677               mr_ozone,cldtau, cldemi)
[3778]1678      ENDIF
1679#endif
[3512]1680
1681#ifdef CPP_COSP2
1682        IF (ok_cosp) THEN
[3726]1683!           DO k = 1, klev
1684!             DO i = 1, klon
1685!               phicosp(i,k) = pphi(i,k) + pphis(i)
1686!             ENDDO
1687!           ENDDO
[3512]1688          CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
1689               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1690               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1691               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1692               JrNt,ref_liq,ref_ice, &
1693               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1694               zu10m,zv10m,pphis, &
1695               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1696               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1697               prfl(:,1:klev),psfl(:,1:klev), &
1698               pmflxr(:,1:klev),pmflxs(:,1:klev), &
1699               mr_ozone,cldtau, cldemi)
1700       ENDIF
1701#endif
1702
1703#ifdef CPP_COSPV2
1704        IF (ok_cosp) THEN
1705           DO k = 1, klev
1706             DO i = 1, klon
1707               phicosp(i,k) = pphi(i,k) + pphis(i)
1708             ENDDO
1709           ENDDO
1710          CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
1711               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1712               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1713               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1714               JrNt,ref_liq,ref_ice, &
1715               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1716               zu10m,zv10m,pphis, &
1717               phicosp,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1718               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1719               prfl(:,1:klev),psfl(:,1:klev), &
1720               pmflxr(:,1:klev),pmflxs(:,1:klev), &
1721               mr_ozone,cldtau, cldemi)
1722       ENDIF
1723#endif
1724
[3460]1725       !
[3465]1726       !
1727!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3460]1728       ! Nouvelle initialisation pour le rayonnement RRTM
[3465]1729!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3435]1730
[3460]1731       CALL iniradia(klon,klev,paprs(1,1:klev+1))
[4085]1732
1733!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1734       CALL wake_ini(rg,rd,rv,prt_level)
[4089]1735       CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
1736   &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
1737
[4085]1738!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1739
[3956]1740       !
1741!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1742       ! Initialisation des champs dans phytrac* qui sont utilises par phys_output_write*
1743       !
1744!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3776]1745
1746#ifdef CPP_Dust
1747       ! Quand on utilise SPLA, on force iflag_phytrac=1
1748       CALL phytracr_spl_out_init()
1749       CALL phys_output_write_spl(itap, pdtphys, paprs, pphis,                  &
1750                                pplay, lmax_th, aerosol_couple,                 &
1751                                ok_ade, ok_aie, ivap, ok_sync,                  &
1752                                ptconv, read_climoz, clevSTD,                   &
1753                                ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse,      &
1754                                flag_aerosol, flag_aerosol_strat, ok_cdnc)
1755#else
1756       ! phys_output_write écrit des variables traceurs seulement si iflag_phytrac == 1
1757       ! donc seulement dans ce cas on doit appeler phytrac_init()
[3465]1758       IF (iflag_phytrac == 1 ) THEN
1759          CALL phytrac_init()
[3776]1760       ENDIF
[3465]1761       CALL phys_output_write(itap, pdtphys, paprs, pphis,                    &
1762                              pplay, lmax_th, aerosol_couple,                 &
[3630]1763                              ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ok_sync,&
[3465]1764                              ptconv, read_climoz, clevSTD,                   &
1765                              ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
1766                              flag_aerosol, flag_aerosol_strat, ok_cdnc)
[3776]1767#endif
[3460]1768
[3776]1769
[3435]1770#ifdef CPP_XIOS
1771       IF (is_omp_master) CALL xios_update_calendar(1)
1772#endif
[3465]1773       IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
[3877]1774       CALL create_etat0_limit_unstruct
1775       CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
[3435]1776
[3465]1777!jyg<
[3577]1778       IF (iflag_pbl<=1) THEN
1779          ! No TKE for Standard Physics
1780          pbl_tke(:,:,:)=0.
1781
1782       ELSE IF (klon_glo==1) THEN
[3465]1783          pbl_tke(:,:,is_ave) = 0.
1784          DO nsrf=1,nbsrf
1785            DO k = 1,klev+1
1786                 pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
1787                     +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
1788            ENDDO
1789          ENDDO
[3988]1790       ELSE
[3465]1791          pbl_tke(:,:,is_ave) = 0. !ym missing init : maybe must be initialized in the same way that for klon_glo==1 ??
1792!>jyg
1793       ENDIF
[2469]1794       !IM begin
1795       print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1796            ,ratqs(1,1)
1797       !IM end
[878]1798
1799
[2469]1800       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1801       !
1802       ! on remet le calendrier a zero
1803       !
1804       IF (raz_date .eq. 1) THEN
1805          itau_phy = 0
1806       ENDIF
[524]1807
[3435]1808!       IF (ABS(phys_tstep-pdtphys).GT.0.001) THEN
1809!          WRITE(lunout,*) 'Pas physique n est pas correct',phys_tstep, &
1810!               pdtphys
1811!          abort_message='Pas physique n est pas correct '
1812!          !           call abort_physic(modname,abort_message,1)
1813!          phys_tstep=pdtphys
1814!       ENDIF
[2469]1815       IF (nlon .NE. klon) THEN
1816          WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,  &
1817               klon
1818          abort_message='nlon et klon ne sont pas coherents'
[2692]1819          CALL abort_physic(modname,abort_message,1)
[2469]1820       ENDIF
1821       IF (nlev .NE. klev) THEN
1822          WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1823               klev
1824          abort_message='nlev et klev ne sont pas coherents'
[2692]1825          CALL abort_physic(modname,abort_message,1)
[2469]1826       ENDIF
1827       !
[3435]1828       IF (phys_tstep*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
[2469]1829          WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1830          WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1831          abort_message='Nbre d appels au rayonnement insuffisant'
[2692]1832          CALL abort_physic(modname,abort_message,1)
[2469]1833       ENDIF
[3956]1834
1835!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1836       ! Initialisation pour la convection de K.E. et pour les poches froides
1837       !
1838!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1839
[2469]1840       WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
[3989]1841       WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", ok_cvl
[2469]1842       !
1843       !KE43
1844       ! Initialisation pour la convection de K.E. (sb):
1845       IF (iflag_con.GE.3) THEN
[524]1846
[2469]1847          WRITE(lunout,*)"*** Convection de Kerry Emanuel 4.3  "
1848          WRITE(lunout,*) &
1849               "On va utiliser le melange convectif des traceurs qui"
1850          WRITE(lunout,*)"est calcule dans convect4.3"
1851          WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
[524]1852
[2469]1853          DO i = 1, klon
1854             ema_cbmf(i) = 0.
1855             ema_pcb(i)  = 0.
1856             ema_pct(i)  = 0.
1857             !          ema_workcbmf(i) = 0.
1858          ENDDO
1859          !IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1860          DO i = 1, klon
1861             ibas_con(i) = 1
1862             itop_con(i) = 1
1863          ENDDO
1864          !IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1865          !================================================================
1866          !CR:04.12.07: initialisations poches froides
1867          ! Controle de ALE et ALP pour la fermeture convective (jyg)
[2692]1868          IF (iflag_wake>=1) THEN
[2469]1869             CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1870                  ,alp_bl_prescr, ale_bl_prescr)
1871             ! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1872             !        print*,'apres ini_wake iflag_cld_th=', iflag_cld_th
[2638]1873             !
1874             ! Initialize tendencies of wake state variables (for some flag values
1875             ! they are not computed).
1876             d_deltat_wk(:,:) = 0.
1877             d_deltaq_wk(:,:) = 0.
1878             d_deltat_wk_gw(:,:) = 0.
1879             d_deltaq_wk_gw(:,:) = 0.
1880             d_deltat_vdf(:,:) = 0.
1881             d_deltaq_vdf(:,:) = 0.
1882             d_deltat_the(:,:) = 0.
1883             d_deltaq_the(:,:) = 0.
1884             d_deltat_ajs_cv(:,:) = 0.
1885             d_deltaq_ajs_cv(:,:) = 0.
1886             d_s_wk(:) = 0.
1887             d_dens_wk(:) = 0.
[3956]1888          ENDIF  !  (iflag_wake>=1)
[973]1889
[2469]1890          !        do i = 1,klon
1891          !           Ale_bl(i)=0.
1892          !           Alp_bl(i)=0.
1893          !        enddo
[1638]1894
[3435]1895       !ELSE
1896       !   ALLOCATE(tabijGCM(0))
1897       !   ALLOCATE(lonGCM(0), latGCM(0))
1898       !   ALLOCATE(iGCM(0), jGCM(0))
[3956]1899       ENDIF  !  (iflag_con.GE.3)
1900       !
[2469]1901       DO i=1,klon
1902          rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1903       ENDDO
[1863]1904
[2469]1905       !34EK
1906       IF (ok_orodr) THEN
[524]1907
[2469]1908          ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1909          ! FH sans doute a enlever de finitivement ou, si on le
1910          ! garde, l'activer justement quand ok_orodr = false.
1911          ! ce rugoro est utilise par la couche limite et fait double emploi
1912          ! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1913          !           DO i=1,klon
1914          !             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1915          !           ENDDO
1916          ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1917          IF (ok_strato) THEN
1918             CALL SUGWD_strato(klon,klev,paprs,pplay)
1919          ELSE
1920             CALL SUGWD(klon,klev,paprs,pplay)
1921          ENDIF
[1863]1922
[2469]1923          DO i=1,klon
1924             zuthe(i)=0.
1925             zvthe(i)=0.
[2692]1926             IF (zstd(i).gt.10.) THEN
[2469]1927                zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1928                zvthe(i)=(1.-zgam(i))*sin(zthe(i))
[2692]1929             ENDIF
[2469]1930          ENDDO
1931       ENDIF
1932       !
1933       !
[3435]1934       lmt_pas = NINT(86400./phys_tstep * 1.0)   ! tous les jours
[2469]1935       WRITE(lunout,*)'La frequence de lecture surface est de ',  &
1936            lmt_pas
1937       !
1938       capemaxcels = 't_max(X)'
1939       t2mincels = 't_min(X)'
1940       t2maxcels = 't_max(X)'
1941       tinst = 'inst(X)'
1942       tave = 'ave(X)'
1943       !IM cf. AM 081204 BEG
[3465]1944       write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
[2469]1945       !IM cf. AM 081204 END
1946       !
1947       !=============================================================
1948       !   Initialisation des sorties
1949       !=============================================================
1950
[3435]1951#ifdef CPP_XIOS
1952       ! Get "missing_val" value from XML files (from temperature variable)
1953       !$OMP MASTER
1954       CALL xios_get_field_attr("temp",default_value=missing_val_omp)
1955       !$OMP END MASTER
1956       !$OMP BARRIER
1957       missing_val=missing_val_omp
1958#endif
1959
[2679]1960#ifdef CPP_XIOS
[3778]1961! Need to put this initialisation after phyetat0 as in the coupled model the XIOS context is only
1962! initialised at that moment
1963       ! Get "missing_val" value from XML files (from temperature variable)
1964       !$OMP MASTER
1965       CALL xios_get_field_attr("temp",default_value=missing_val_omp)
1966       !$OMP END MASTER
1967       !$OMP BARRIER
1968       missing_val=missing_val_omp
[3988]1969       !
1970       ! Now we activate some double radiation call flags only if some
1971       ! diagnostics are requested, otherwise there is no point in doing this
1972       IF (is_master) THEN
1973         !--setting up swaero_diag to TRUE in XIOS case
1974         IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
1975            xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
1976            xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR.  &
1977              (iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
1978                                  xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0"))))  &
1979            !!!--for now these fields are not in the XML files so they are omitted
1980            !!!  xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
1981            swaero_diag=.TRUE.
1982 
1983         !--setting up swaerofree_diag to TRUE in XIOS case
1984         IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
1985            xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR.   &
1986            xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
1987            xios_field_is_active("LWupTOAcleanclr")) &
1988            swaerofree_diag=.TRUE.
1989 
1990         !--setting up dryaod_diag to TRUE in XIOS case
1991         DO naero = 1, naero_tot-1
1992          IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
1993         ENDDO
1994         !
1995         !--setting up ok_4xCO2atm to TRUE in XIOS case
1996         IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
1997            xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
1998            xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
1999            xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
2000            xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
2001            xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
2002            ok_4xCO2atm=.TRUE.
2003       ENDIF
2004       !$OMP BARRIER
2005       CALL bcast(swaero_diag)
2006       CALL bcast(swaerofree_diag)
2007       CALL bcast(dryaod_diag)
2008       CALL bcast(ok_4xCO2atm)
[2679]2009#endif
[3988]2010       !
[3435]2011       CALL printflag( tabcntr0,radpas,ok_journe, &
2012            ok_instan, ok_region )
[2469]2013       !
2014       !
2015       ! Prescrire l'ozone dans l'atmosphere
2016       !
2017       !c         DO i = 1, klon
2018       !c         DO k = 1, klev
2019       !c            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
2020       !c         ENDDO
2021       !c         ENDDO
2022       !
[3870]2023       IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN                   ! ModThL
[524]2024#ifdef INCA
[2469]2025          CALL VTe(VTphysiq)
2026          CALL VTb(VTinca)
2027          calday = REAL(days_elapsed) + jH_cur
2028          WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
[959]2029
[2469]2030          CALL chemini(  &
2031               rg, &
2032               ra, &
2033               cell_area, &
2034               latitude_deg, &
2035               longitude_deg, &
2036               presnivs, &
2037               calday, &
2038               klon, &
2039               nqtot, &
[3872]2040               nqo+nqCO2, &
[2469]2041               pdtphys, &
2042               annee_ref, &
[2906]2043               year_cur, &
[2469]2044               day_ref,  &
2045               day_ini, &
2046               start_time, &
2047               itau_phy, &
[2831]2048               date0, &
[2469]2049               io_lon, &
[3351]2050               io_lat, &
[3418]2051               chemistry_couple, &
2052               init_source, &
2053               init_tauinca, &
2054               init_pizinca, &
2055               init_cginca, &
2056               init_ccminca)
[959]2057
[3418]2058
2059          ! initialisation des variables depuis le restart de inca
2060          ccm(:,:,:) = init_ccminca
2061          tau_aero(:,:,:,:) = init_tauinca
2062          piz_aero(:,:,:,:) = init_pizinca
2063          cg_aero(:,:,:,:) = init_cginca
2064!         
2065
2066
[2469]2067          CALL VTe(VTinca)
2068          CALL VTb(VTphysiq)
[524]2069#endif
[2692]2070       ENDIF
[3988]2071       !
[3666]2072       IF (type_trac == 'repr') THEN
2073#ifdef REPROBUS
2074          CALL chemini_rep(  &
2075               presnivs, &
2076               pdtphys, &
2077               annee_ref, &
2078               day_ref,  &
2079               day_ini, &
2080               start_time, &
2081               itau_phy, &
2082               io_lon, &
2083               io_lat)
2084#endif
2085       ENDIF
[3465]2086
[2469]2087       !$omp single
[2788]2088       IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz,   &
[2820]2089           press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
[2469]2090       !$omp end single
2091       !
2092       !IM betaCRF
2093       pfree=70000. !Pa
2094       beta_pbl=1.
2095       beta_free=1.
2096       lon1_beta=-180.
2097       lon2_beta=+180.
2098       lat1_beta=90.
2099       lat2_beta=-90.
2100       mskocean_beta=.FALSE.
[1539]2101
[2469]2102       !albedo SB >>>
[3461]2103       SELECT CASE(nsw)
2104       CASE(2)
[2469]2105          SFRWL(1)=0.45538747
2106          SFRWL(2)=0.54461211
[3461]2107       CASE(4)
[2469]2108          SFRWL(1)=0.45538747
2109          SFRWL(2)=0.32870591
2110          SFRWL(3)=0.18568763
2111          SFRWL(4)=3.02191470E-02
[3461]2112       CASE(6)
[2469]2113          SFRWL(1)=1.28432794E-03
2114          SFRWL(2)=0.12304168
2115          SFRWL(3)=0.33106142
2116          SFRWL(4)=0.32870591
2117          SFRWL(5)=0.18568763
2118          SFRWL(6)=3.02191470E-02
[3461]2119       END SELECT
[2469]2120       !albedo SB <<<
[2227]2121
[2469]2122       OPEN(99,file='beta_crf.data',status='old', &
2123            form='formatted',err=9999)
2124       READ(99,*,end=9998) pfree
2125       READ(99,*,end=9998) beta_pbl
2126       READ(99,*,end=9998) beta_free
2127       READ(99,*,end=9998) lon1_beta
2128       READ(99,*,end=9998) lon2_beta
2129       READ(99,*,end=9998) lat1_beta
2130       READ(99,*,end=9998) lat2_beta
2131       READ(99,*,end=9998) mskocean_beta
21329998   Continue
2133       CLOSE(99)
21349999   Continue
2135       WRITE(*,*)'pfree=',pfree
2136       WRITE(*,*)'beta_pbl=',beta_pbl
2137       WRITE(*,*)'beta_free=',beta_free
2138       WRITE(*,*)'lon1_beta=',lon1_beta
2139       WRITE(*,*)'lon2_beta=',lon2_beta
2140       WRITE(*,*)'lat1_beta=',lat1_beta
2141       WRITE(*,*)'lat2_beta=',lat2_beta
2142       WRITE(*,*)'mskocean_beta=',mskocean_beta
[3048]2143
2144      !lwoff=y : offset LW CRE for radiation code and other schemes
2145      !lwoff=y : betalwoff=1.
2146      betalwoff=0.
2147      IF (ok_lwoff) THEN
2148         betalwoff=1.
2149      ENDIF
2150      WRITE(*,*)'ok_lwoff=',ok_lwoff
2151      !
2152      !lwoff=y to begin only sollw and sollwdown are set up to CS values
2153      sollw = sollw + betalwoff * (sollw0 - sollw)
2154      sollwdown(:)= sollwdown(:) + betalwoff *(-1.*ZFLDN0(:,1) - &
2155                    sollwdown(:))
[3597]2156
2157
[4085]2158
[2469]2159    ENDIF
2160    !
2161    !   ****************     Fin  de   IF ( debut  )   ***************
2162    !
2163    !
2164    ! Incrementer le compteur de la physique
2165    !
2166    itap   = itap + 1
[2795]2167    IF (is_master .OR. prt_level > 9) THEN
[2783]2168      IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
[2795]2169         WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
2170         WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
2171 100     FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
[2783]2172      ENDIF
2173    ENDIF
[2469]2174    !
2175    !
2176    ! Update fraction of the sub-surfaces (pctsrf) and
2177    ! initialize, where a new fraction has appeared, all variables depending
2178    ! on the surface fraction.
2179    !
[3435]2180    CALL change_srf_frac(itap, phys_tstep, days_elapsed+1,  &
[2469]2181         pctsrf, fevap, z0m, z0h, agesno,              &
2182         falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
[996]2183
[2469]2184    ! Update time and other variables in Reprobus
2185    IF (type_trac == 'repr') THEN
[1565]2186#ifdef REPROBUS
[2469]2187       CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
2188       print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
2189       CALL Rtime(debut)
[1565]2190#endif
[2692]2191    ENDIF
[1565]2192
[2469]2193    ! Tendances bidons pour les processus qui n'affectent pas certaines
2194    ! variables.
2195    du0(:,:)=0.
2196    dv0(:,:)=0.
2197    dt0 = 0.
2198    dq0(:,:)=0.
2199    dql0(:,:)=0.
2200    dqi0(:,:)=0.
[2635]2201    dsig0(:) = 0.
2202    ddens0(:) = 0.
2203    wkoccur1(:)=1
[2469]2204    !
2205    ! Mettre a zero des variables de sortie (pour securite)
2206    !
2207    DO i = 1, klon
2208       d_ps(i) = 0.0
2209    ENDDO
2210    DO k = 1, klev
2211       DO i = 1, klon
2212          d_t(i,k) = 0.0
2213          d_u(i,k) = 0.0
2214          d_v(i,k) = 0.0
2215       ENDDO
2216    ENDDO
2217    DO iq = 1, nqtot
2218       DO k = 1, klev
2219          DO i = 1, klon
2220             d_qx(i,k,iq) = 0.0
2221          ENDDO
2222       ENDDO
2223    ENDDO
2224    beta_prec_fisrt(:,:)=0.
2225    beta_prec(:,:)=0.
[3134]2226    !
2227    !   Output variables from the convective scheme should not be set to 0
2228    !   since convection is not always called at every time step.
2229    IF (ok_bug_cv_trac) THEN
2230      da(:,:)=0.
2231      mp(:,:)=0.
2232      phi(:,:,:)=0.
2233      ! RomP >>>
2234      phi2(:,:,:)=0.
2235      epmlmMm(:,:,:)=0.
2236      eplaMm(:,:)=0.
2237      d1a(:,:)=0.
2238      dam(:,:)=0.
2239      pmflxr(:,:)=0.
2240      pmflxs(:,:)=0.
2241      ! RomP <<<
2242    ENDIF
[2469]2243    !
2244    ! Ne pas affecter les valeurs entrees de u, v, h, et q
2245    !
2246    DO k = 1, klev
2247       DO i = 1, klon
2248          t_seri(i,k)  = t(i,k)
2249          u_seri(i,k)  = u(i,k)
2250          v_seri(i,k)  = v(i,k)
2251          q_seri(i,k)  = qx(i,k,ivap)
2252          ql_seri(i,k) = qx(i,k,iliq)
2253          !CR: ATTENTION, on rajoute la variable glace
[4098]2254          IF (nqo.EQ.2) THEN             !--vapour and liquid only
[2469]2255             qs_seri(i,k) = 0.
[4066]2256             rneb_seri(i,k) = 0.
[4098]2257          ELSE IF (nqo.EQ.3) THEN        !--vapour, liquid and ice
[2469]2258             qs_seri(i,k) = qx(i,k,isol)
[4066]2259             rneb_seri(i,k) = 0.
[4098]2260          ELSE IF (nqo.EQ.4) THEN        !--vapour, liquid, ice and rneb
[4059]2261             qs_seri(i,k) = qx(i,k,isol)
2262             rneb_seri(i,k) = qx(i,k,irneb)
[2692]2263          ENDIF
[2469]2264       ENDDO
2265    ENDDO
[2476]2266    !
2267    !--OB mass fixer
2268    IF (mass_fixer) THEN
2269    !--store initial water burden
2270    qql1(:)=0.0
[2499]2271    DO k = 1, klev
2272      qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
[2476]2273    ENDDO
2274    ENDIF
2275    !--fin mass fixer
2276
[2469]2277    tke0(:,:)=pbl_tke(:,:,is_ave)
[4056]2278    IF (nqtot > nqo) THEN
2279       ! water isotopes are not included in tr_seri
2280       itr = 0
2281       DO iq = 1, nqtot
[4071]2282         IF(.NOT.tracers(iq)%isInPhysics) CYCLE
[4056]2283         itr = itr+1
[2469]2284          DO  k = 1, klev
2285             DO  i = 1, klon
[4056]2286                tr_seri(i,k,itr) = qx(i,k,iq)
[2469]2287             ENDDO
2288          ENDDO
2289       ENDDO
2290    ELSE
[4056]2291! DC: make sure the final "1" index was meant for 1st H2O phase (vapor) !!!
[4069]2292!       tr_seri(:,:,strIdx(tracers(:)%name,addPhase('H2O','g'))) = 0.0
2293       tr_seri(:,:,strIdx(tracers(:)%name,addPhase('H2O','v',''))) = 0.0
[2469]2294    ENDIF
[3599]2295!
2296! Temporary solutions adressing ticket #104 and the non initialisation of tr_ancien
2297! LF
2298    IF (debut) THEN
2299      WRITE(lunout,*)' WARNING: tr_ancien initialised to tr_seri'
[4056]2300       itr = 0
2301       do iq = 1, nqtot
[4071]2302         IF(.NOT.tracers(iq)%isInPhysics) CYCLE
[4056]2303         itr = itr+1
2304         tr_ancien(:,:,itr)=tr_seri(:,:,itr)       
2305       enddo
[3599]2306    ENDIF
[2469]2307    !
2308    DO i = 1, klon
2309       ztsol(i) = 0.
2310    ENDDO
2311    DO nsrf = 1, nbsrf
2312       DO i = 1, klon
2313          ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
2314       ENDDO
2315    ENDDO
[2611]2316    ! Initialize variables used for diagnostic purpose
[2692]2317    IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
[524]2318
[2469]2319    ! Diagnostiquer la tendance dynamique
2320    !
2321    IF (ancien_ok) THEN
[2499]2322    !
[3435]2323       d_u_dyn(:,:)  = (u_seri(:,:)-u_ancien(:,:))/phys_tstep
2324       d_v_dyn(:,:)  = (v_seri(:,:)-v_ancien(:,:))/phys_tstep
2325       d_t_dyn(:,:)  = (t_seri(:,:)-t_ancien(:,:))/phys_tstep
2326       d_q_dyn(:,:)  = (q_seri(:,:)-q_ancien(:,:))/phys_tstep
2327       d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
2328       d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
[2499]2329       CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
[3435]2330       d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
[2499]2331       CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
[3435]2332       d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/phys_tstep
[2499]2333       CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
[3435]2334       d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
[2469]2335       ! !! RomP >>>   td dyn traceur
[4056]2336       IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
[2469]2337       ! !! RomP <<<
[4059]2338       !!d_rneb_dyn(:,:)=(rneb_seri(:,:)-rneb_ancien(:,:))/phys_tstep
2339       d_rneb_dyn(:,:)=0.0
[2469]2340    ELSE
[2499]2341       d_u_dyn(:,:)  = 0.0
2342       d_v_dyn(:,:)  = 0.0
2343       d_t_dyn(:,:)  = 0.0
2344       d_q_dyn(:,:)  = 0.0
2345       d_ql_dyn(:,:) = 0.0
2346       d_qs_dyn(:,:) = 0.0
2347       d_q_dyn2d(:)  = 0.0
2348       d_ql_dyn2d(:) = 0.0
2349       d_qs_dyn2d(:) = 0.0
[2469]2350       ! !! RomP >>>   td dyn traceur
[4056]2351       IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
[2469]2352       ! !! RomP <<<
[4059]2353       d_rneb_dyn(:,:)=0.0
[2469]2354       ancien_ok = .TRUE.
2355    ENDIF
2356    !
2357    ! Ajouter le geopotentiel du sol:
2358    !
2359    DO k = 1, klev
2360       DO i = 1, klon
2361          zphi(i,k) = pphi(i,k) + pphis(i)
2362       ENDDO
2363    ENDDO
2364    !
2365    ! Verifier les temperatures
2366    !
2367    !IM BEG
2368    IF (check) THEN
2369       amn=MIN(ftsol(1,is_ter),1000.)
2370       amx=MAX(ftsol(1,is_ter),-1000.)
2371       DO i=2, klon
2372          amn=MIN(ftsol(i,is_ter),amn)
2373          amx=MAX(ftsol(i,is_ter),amx)
2374       ENDDO
2375       !
2376       PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
2377    ENDIF !(check) THEN
2378    !IM END
2379    !
2380    CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
2381    IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
[2235]2382
[2469]2383    !
2384    !IM BEG
2385    IF (check) THEN
2386       amn=MIN(ftsol(1,is_ter),1000.)
2387       amx=MAX(ftsol(1,is_ter),-1000.)
2388       DO i=2, klon
2389          amn=MIN(ftsol(i,is_ter),amn)
2390          amx=MAX(ftsol(i,is_ter),amx)
2391       ENDDO
2392       !
2393       PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2394    ENDIF !(check) THEN
2395    !IM END
2396    !
2397    ! Mettre en action les conditions aux limites (albedo, sst, etc.).
2398    ! Prescrire l'ozone et calculer l'albedo sur l'ocean.
2399    !
[2661]2400    ! Update ozone if day change
2401    IF (MOD(itap-1,lmt_pas) == 0) THEN
[2774]2402       IF (read_climoz <= 0) THEN
2403          ! Once per day, update ozone from Royer:
2404          IF (solarlong0<-999.) then
2405             ! Generic case with evolvoing season
2406             zzz=real(days_elapsed+1)
2407          ELSE IF (abs(solarlong0-1000.)<1.e-4) then
2408             ! Particular case with annual mean insolation
2409             zzz=real(90) ! could be revisited
2410             IF (read_climoz/=-1) THEN
2411                abort_message ='read_climoz=-1 is recommended when ' &
2412                     // 'solarlong0=1000.'
2413                CALL abort_physic (modname,abort_message,1)
2414             ENDIF
2415          ELSE
2416             ! Case where the season is imposed with solarlong0
2417             zzz=real(90) ! could be revisited
2418          ENDIF
[2661]2419
[2774]2420          wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
[3666]2421#ifdef REPROBUS
2422          ptrop=dyn_tropopause(t_seri, ztsol, paprs, pplay, rot)/100.
2423          DO i = 1, klon
2424             Z1=t_seri(i,itroprep(i)+1)
2425             Z2=t_seri(i,itroprep(i))
2426             fac=(Z1-Z2)/alog(pplay(i,itroprep(i)+1)/pplay(i,itroprep(i)))
2427             B=Z2-fac*alog(pplay(i,itroprep(i)))
2428             ttrop(i)= fac*alog(ptrop(i))+B
2429!       
2430             Z1= 1.e-3 * ( pphi(i,itroprep(i)+1)+pphis(i) ) / gravit
2431             Z2= 1.e-3 * ( pphi(i,itroprep(i))  +pphis(i) ) / gravit
2432             fac=(Z1-Z2)/alog(pplay(i,itroprep(i)+1)/pplay(i,itroprep(i)))
2433             B=Z2-fac*alog(pplay(i,itroprep(i)))
2434             ztrop(i)=fac*alog(ptrop(i))+B
2435          ENDDO
2436#endif
[2774]2437       ELSE
[2820]2438          !--- ro3i = elapsed days number since current year 1st january, 0h
2439          ro3i=days_elapsed+jh_cur-jh_1jan
2440          !--- scaling for old style files (360 records)
2441          IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
[2788]2442          IF(adjust_tropopause) THEN
[2820]2443             CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),   &
[3086]2444                      ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1),    &
2445                      time_climoz ,  longitude_deg,   latitude_deg,          &
[2968]2446                      dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
[2774]2447          ELSE
[3086]2448             CALL regr_pr_time_av(ncid_climoz,  vars_climoz(1:read_climoz),  &
2449                      ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1),    &
2450                      time_climoz )
[3461]2451          ENDIF
[2774]2452          ! Convert from mole fraction of ozone to column density of ozone in a
2453          ! cell, in kDU:
2454          FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
2455               * zmasse / dobson_u / 1e3
[2788]2456          ! (By regridding ozone values for LMDZ only once a day, we
[2774]2457          ! have already neglected the variation of pressure in one
2458          ! day. So do not recompute "wo" at each time step even if
2459          ! "zmasse" changes a little.)
2460       ENDIF
[2469]2461    ENDIF
2462    !
2463    ! Re-evaporer l'eau liquide nuageuse
2464    !
[2705]2465     CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
2466   &         d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
[1849]2467
[2705]2468     CALL add_phys_tend &
2469            (du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
[2812]2470               'eva',abortphy,flag_inhib_tend,itap,0)
[3461]2471    CALL prt_enerbil('eva',itap)
[2086]2472
[2469]2473    !=========================================================================
2474    ! Calculs de l'orbite.
2475    ! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2476    ! doit donc etre plac\'e avant radlwsw et pbl_surface
[883]2477
[2469]2478    ! !!   jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2692]2479    CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
[2469]2480    day_since_equinox = (jD_cur + jH_cur) - jD_eq
2481    !
2482    !   choix entre calcul de la longitude solaire vraie ou valeur fixee a
2483    !   solarlong0
[2692]2484    IF (solarlong0<-999.) THEN
2485       IF (new_orbit) THEN
[2469]2486          ! calcul selon la routine utilisee pour les planetes
[2692]2487          CALL solarlong(day_since_equinox, zlongi, dist)
2488       ELSE
[2469]2489          ! calcul selon la routine utilisee pour l'AR4
2490          CALL orbite(REAL(days_elapsed+1),zlongi,dist)
[2692]2491       ENDIF
2492    ELSE
[2469]2493       zlongi=solarlong0  ! longitude solaire vraie
2494       dist=1.            ! distance au soleil / moyenne
[2692]2495    ENDIF
[1529]2496
[2692]2497    IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
[524]2498
[2692]2499
[2469]2500    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2501    ! Calcul de l'ensoleillement :
2502    ! ============================
2503    ! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2504    ! l'annee a partir d'une formule analytique.
2505    ! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2506    ! non nul aux poles.
[2692]2507    IF (abs(solarlong0-1000.)<1.e-4) THEN
2508       CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
[2469]2509            latitude_deg,longitude_deg,rmu0,fract)
[2979]2510       swradcorr(:) = 1.0
2511       JrNt(:) = 1.0
2512       zrmu0(:) = rmu0(:)
[2469]2513    ELSE
2514       ! recode par Olivier Boucher en sept 2015
2515       SELECT CASE (iflag_cycle_diurne)
2516       CASE(0) 
2517          !  Sans cycle diurne
2518          CALL angle(zlongi, latitude_deg, fract, rmu0)
2519          swradcorr = 1.0
2520          JrNt = 1.0
2521          zrmu0 = rmu0
2522       CASE(1) 
2523          !  Avec cycle diurne sans application des poids
2524          !  bit comparable a l ancienne formulation cycle_diurne=true
2525          !  on integre entre gmtime et gmtime+radpas
[3435]2526          zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s)
[2469]2527          CALL zenang(zlongi,jH_cur,0.0,zdtime, &
2528               latitude_deg,longitude_deg,rmu0,fract)
2529          zrmu0 = rmu0
2530          swradcorr = 1.0
2531          ! Calcul du flag jour-nuit
2532          JrNt = 0.0
2533          WHERE (fract.GT.0.0) JrNt = 1.0
2534       CASE(2) 
2535          !  Avec cycle diurne sans application des poids
2536          !  On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
2537          !  Comme cette routine est appele a tous les pas de temps de
2538          !  la physique meme si le rayonnement n'est pas appele je
2539          !  remonte en arriere les radpas-1 pas de temps
2540          !  suivant. Petite ruse avec MOD pour prendre en compte le
2541          !  premier pas de temps de la physique pendant lequel
2542          !  itaprad=0
[3435]2543          zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1)     
2544          zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1)
[2469]2545          CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2546               latitude_deg,longitude_deg,rmu0,fract)
2547          !
2548          ! Calcul des poids
2549          !
[3435]2550          zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le
[2469]2551          zdtime2=0.0    !--pas de temps de la physique qui se termine
2552          CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2553               latitude_deg,longitude_deg,zrmu0,zfract)
2554          swradcorr = 0.0
2555          WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
2556               swradcorr=zfract/fract*zrmu0/rmu0
2557          ! Calcul du flag jour-nuit
2558          JrNt = 0.0
2559          WHERE (zfract.GT.0.0) JrNt = 1.0
2560       END SELECT
2561    ENDIF
[3110]2562    sza_o = ACOS (rmu0) *180./pi
[782]2563
[2692]2564    IF (mydebug) THEN
2565       CALL writefield_phy('u_seri',u_seri,nbp_lev)
2566       CALL writefield_phy('v_seri',v_seri,nbp_lev)
2567       CALL writefield_phy('t_seri',t_seri,nbp_lev)
2568       CALL writefield_phy('q_seri',q_seri,nbp_lev)
2569    ENDIF
[883]2570
[2469]2571    !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2572    ! Appel au pbl_surface : Planetary Boudary Layer et Surface
2573    ! Cela implique tous les interactions des sous-surfaces et la
2574    ! partie diffusion turbulent du couche limit.
2575    !
2576    ! Certains varibales de sorties de pbl_surface sont utiliser que pour
2577    ! ecriture des fihiers hist_XXXX.nc, ces sont :
2578    !   qsol,      zq2m,      s_pblh,  s_lcl,
2579    !   s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
2580    !   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
2581    !   zu10m,     zv10m,   fder,
[3888]2582    !   zxqsurf,   delta_qsurf,
2583    !   rh2m,      zxfluxu, zxfluxv,
[2469]2584    !   frugs,     agesno,    fsollw,  fsolsw,
2585    !   d_ts,      fevap,     fluxlat, t2m,
2586    !   wfbils,    wfbilo,    fluxt,   fluxu, fluxv,
2587    !
2588    ! Certains ne sont pas utiliser du tout :
2589    !   dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2590    !
[1724]2591
[2469]2592    ! Calcul de l'humidite de saturation au niveau du sol
[1724]2593
2594
[996]2595
[2692]2596    IF (iflag_pbl/=0) THEN
[2240]2597
[2469]2598       !jyg+nrlmd<
[2852]2599!!jyg       IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
2600       IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,10) .ge. 1) THEN
[3179]2601          print *,'debut du splitting de la PBL, wake_s = ', wake_s(:)
2602          print *,'debut du splitting de la PBL, wake_deltat = ', wake_deltat(:,1)
2603          print *,'debut du splitting de la PBL, wake_deltaq = ', wake_deltaq(:,1)
[2469]2604       ENDIF
2605       ! !!
2606       !>jyg+nrlmd
2607       !
2608       !-------gustiness calculation-------!
[3435]2609       !ym : Warning gustiness non inialized for iflag_gusts=2 & iflag_gusts=3
2610       gustiness=0  !ym missing init
2611       
[2469]2612       IF (iflag_gusts==0) THEN
2613          gustiness(1:klon)=0
2614       ELSE IF (iflag_gusts==1) THEN
[3111]2615          gustiness(1:klon)=f_gust_bl*ale_bl(1:klon)+f_gust_wk*ale_wake(1:klon)
2616       ELSE IF (iflag_gusts==2) THEN
2617          gustiness(1:klon)=f_gust_bl*ale_bl_stat(1:klon)+f_gust_wk*ale_wake(1:klon)
[2469]2618          ! ELSE IF (iflag_gusts==2) THEN
2619          !    do i = 1, klon
2620          !       gustiness(i)=f_gust_bl*ale_bl(i)+sigma_wk(i)*f_gust_wk&
2621          !           *ale_wake(i) !! need to make sigma_wk accessible here
2622          !    enddo
2623          ! ELSE IF (iflag_gusts==3) THEN
2624          !    do i = 1, klon
2625          !       gustiness(i)=f_gust_bl*alp_bl(i)+f_gust_wk*alp_wake(i)
2626          !    enddo
2627       ENDIF
[2278]2628
[2469]2629       CALL pbl_surface(  &
[3435]2630            phys_tstep,     date0,     itap,    days_elapsed+1, &
[2469]2631            debut,     lafin, &
2632            longitude_deg, latitude_deg, rugoro,  zrmu0,      &
[4056]2633            sollwdown,    cldt,      &
[3756]2634            rain_fall, snow_fall, solsw,   solswfdiff, sollw,     &
[2469]2635            gustiness,                                &
2636            t_seri,    q_seri,    u_seri,  v_seri,    &
2637                                !nrlmd+jyg<
2638            wake_deltat, wake_deltaq, wake_cstar, wake_s, &
2639                                !>nrlmd+jyg
2640            pplay,     paprs,     pctsrf,             &
2641            ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
2642                                !albedo SB <<<
2643            cdragh,    cdragm,  u1,    v1,            &
[3888]2644            beta_aridity, &
[2469]2645                                !albedo SB >>>
2646                                ! albsol1,   albsol2,   sens,    evap,      &
2647            albsol_dir,   albsol_dif,   sens,    evap,   & 
2648                                !albedo SB <<<
2649            albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
[3817]2650            zxtsol,    zxfluxlat, zt2m,    qsat2m,  zn2mout, &
[2469]2651            d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss, &
2652                                !nrlmd<
2653                                !jyg<
2654            d_t_vdf_w, d_q_vdf_w, &
2655            d_t_vdf_x, d_q_vdf_x, &
2656            sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
2657                                !>jyg
2658            delta_tsurf,wake_dens, &
2659            cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
2660            kh,kh_x,kh_w, &
2661                                !>nrlmd
2662            coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
2663            slab_wfbils,                 &
2664            qsol,      zq2m,      s_pblh,  s_lcl, &
2665                                !jyg<
2666            s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
2667                                !>jyg
2668            s_capCL,   s_oliqCL,  s_cteiCL,s_pblT, &
2669            s_therm,   s_trmb1,   s_trmb2, s_trmb3, &
2670            zustar, zu10m,     zv10m,   fder, &
[3888]2671            zxqsurf, delta_qsurf,   rh2m,      zxfluxu, zxfluxv, &
[2469]2672            z0m, z0h,     agesno,    fsollw,  fsolsw, &
2673            d_ts,      fevap,     fluxlat, t2m, &
[2670]2674            wfbils, wfbilo, wfevap, wfrain, wfsnow, &
2675            fluxt,   fluxu,  fluxv, &
[2469]2676            dsens,     devap,     zxsnow, &
2677            zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke, &
2678                                !nrlmd+jyg<
[2952]2679            wake_delta_pbl_TKE, &
[2469]2680                                !>nrlmd+jyg
[2952]2681             treedrg )
2682!FC
[2469]2683       !
2684       !  Add turbulent diffusion tendency to the wake difference variables
[2852]2685!!jyg       IF (mod(iflag_pbl_split,2) .NE. 0) THEN
2686       IF (mod(iflag_pbl_split,10) .NE. 0) THEN
[2635]2687!jyg<
2688          d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
2689          d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
2690          CALL add_wake_tend &
[3208]2691             (d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy)
[2635]2692       ELSE
2693          d_deltat_vdf(:,:) = 0.
2694          d_deltaq_vdf(:,:) = 0.
2695!>jyg
[2469]2696       ENDIF
[1624]2697
[2469]2698       !---------------------------------------------------------------------
2699       ! ajout des tendances de la diffusion turbulente
2700       IF (klon_glo==1) THEN
2701          CALL add_pbl_tend &
2702               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
[2799]2703               'vdf',abortphy,flag_inhib_tend,itap)
[2469]2704       ELSE
2705          CALL add_phys_tend &
2706               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
[2812]2707               'vdf',abortphy,flag_inhib_tend,itap,0)
[2469]2708       ENDIF
[3461]2709       CALL prt_enerbil('vdf',itap)
[2469]2710       !--------------------------------------------------------------------
[766]2711
[2692]2712       IF (mydebug) THEN
2713          CALL writefield_phy('u_seri',u_seri,nbp_lev)
2714          CALL writefield_phy('v_seri',v_seri,nbp_lev)
2715          CALL writefield_phy('t_seri',t_seri,nbp_lev)
2716          CALL writefield_phy('q_seri',q_seri,nbp_lev)
2717       ENDIF
[2227]2718
[2469]2719       !albedo SB >>>
2720       albsol1=0.
2721       albsol2=0.
2722       falb1=0.
2723       falb2=0.
[2692]2724       SELECT CASE(nsw)
2725       CASE(2)
[2469]2726          albsol1=albsol_dir(:,1)
2727          albsol2=albsol_dir(:,2)
2728          falb1=falb_dir(:,1,:)
2729          falb2=falb_dir(:,2,:)
[2692]2730       CASE(4)
[2469]2731          albsol1=albsol_dir(:,1)
2732          albsol2=albsol_dir(:,2)*SFRWL(2)+albsol_dir(:,3)*SFRWL(3) &
2733               +albsol_dir(:,4)*SFRWL(4)
2734          albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2735          falb1=falb_dir(:,1,:)
2736          falb2=falb_dir(:,2,:)*SFRWL(2)+falb_dir(:,3,:)*SFRWL(3) &
2737               +falb_dir(:,4,:)*SFRWL(4)
2738          falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
[2692]2739       CASE(6)
[2469]2740          albsol1=albsol_dir(:,1)*SFRWL(1)+albsol_dir(:,2)*SFRWL(2) &
2741               +albsol_dir(:,3)*SFRWL(3)
2742          albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2743          albsol2=albsol_dir(:,4)*SFRWL(4)+albsol_dir(:,5)*SFRWL(5) &
2744               +albsol_dir(:,6)*SFRWL(6)
2745          albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2746          falb1=falb_dir(:,1,:)*SFRWL(1)+falb_dir(:,2,:)*SFRWL(2) &
2747               +falb_dir(:,3,:)*SFRWL(3)
2748          falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2749          falb2=falb_dir(:,4,:)*SFRWL(4)+falb_dir(:,5,:)*SFRWL(5) &
2750               +falb_dir(:,6,:)*SFRWL(6)
2751          falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
[2692]2752       END SELECt
[2469]2753       !albedo SB <<<
[2227]2754
[766]2755
[2469]2756       CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
2757            t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
[1724]2758
[2469]2759    ENDIF
2760    ! =================================================================== c
2761    !   Calcul de Qsat
[881]2762
[2469]2763    DO k = 1, klev
2764       DO i = 1, klon
2765          zx_t = t_seri(i,k)
2766          IF (thermcep) THEN
2767             zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2768             zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2769             zx_qs  = MIN(0.5,zx_qs)
2770             zcor   = 1./(1.-retv*zx_qs)
2771             zx_qs  = zx_qs*zcor
2772          ELSE
2773             !!           IF (zx_t.LT.t_coup) THEN             !jyg
2774             IF (zx_t.LT.rtt) THEN                  !jyg
2775                zx_qs = qsats(zx_t)/pplay(i,k)
2776             ELSE
2777                zx_qs = qsatl(zx_t)/pplay(i,k)
2778             ENDIF
2779          ENDIF
2780          zqsat(i,k)=zx_qs
2781       ENDDO
2782    ENDDO
[959]2783
[2692]2784    IF (prt_level.ge.1) THEN
[2469]2785       write(lunout,*) 'L   qsat (g/kg) avant clouds_gno'
2786       write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
[2692]2787    ENDIF
[2469]2788    !
2789    ! Appeler la convection (au choix)
2790    !
2791    DO k = 1, klev
2792       DO i = 1, klon
2793          conv_q(i,k) = d_q_dyn(i,k)  &
[3435]2794               + d_q_vdf(i,k)/phys_tstep
[2469]2795          conv_t(i,k) = d_t_dyn(i,k)  &
[3435]2796               + d_t_vdf(i,k)/phys_tstep
[2469]2797       ENDDO
2798    ENDDO
2799    IF (check) THEN
2800       za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2801       WRITE(lunout,*) "avantcon=", za
2802    ENDIF
2803    zx_ajustq = .FALSE.
2804    IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2805    IF (zx_ajustq) THEN
2806       DO i = 1, klon
2807          z_avant(i) = 0.0
2808       ENDDO
2809       DO k = 1, klev
2810          DO i = 1, klon
2811             z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
2812                  *(paprs(i,k)-paprs(i,k+1))/RG
2813          ENDDO
2814       ENDDO
2815    ENDIF
[959]2816
[2469]2817    ! Calcule de vitesse verticale a partir de flux de masse verticale
2818    DO k = 1, klev
2819       DO i = 1, klon
2820          omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
[2692]2821       ENDDO
2822    ENDDO
2823
2824    IF (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
[2469]2825         omega(igout, :)
[2707]2826    !
2827    ! Appel de la convection tous les "cvpas"
2828    !
[3150]2829!!jyg    IF (MOD(itapcv,cvpas).EQ.0) THEN
[3151]2830!!    print *,' physiq : itapcv, cvpas, itap-1, cvpas_0 ', &
2831!!                       itapcv, cvpas, itap-1, cvpas_0
[3150]2832    IF (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itap-1,cvpas_0).EQ.0) THEN
[2707]2833
[3134]2834    !
2835    ! Mettre a zero des variables de sortie (pour securite)
2836    !
2837    pmflxr(:,:) = 0.
2838    pmflxs(:,:) = 0.
2839    wdtrainA(:,:) = 0.
[3496]2840    wdtrainS(:,:) = 0.
[3134]2841    wdtrainM(:,:) = 0.
2842    upwd(:,:) = 0.
2843    dnwd(:,:) = 0.
2844    ep(:,:) = 0.
2845    da(:,:)=0.
2846    mp(:,:)=0.
2847    wght_cvfd(:,:)=0.
2848    phi(:,:,:)=0.
2849    phi2(:,:,:)=0.
2850    epmlmMm(:,:,:)=0.
2851    eplaMm(:,:)=0.
2852    d1a(:,:)=0.
2853    dam(:,:)=0.
2854    elij(:,:,:)=0.
2855    ev(:,:)=0.
[3496]2856    qtaa(:,:)=0.
[3134]2857    clw(:,:)=0.
2858    sij(:,:,:)=0.
2859    !
[2469]2860    IF (iflag_con.EQ.1) THEN
2861       abort_message ='reactiver le call conlmd dans physiq.F'
2862       CALL abort_physic (modname,abort_message,1)
[3435]2863       !     CALL conlmd (phys_tstep, paprs, pplay, t_seri, q_seri, conv_q,
[2469]2864       !    .             d_t_con, d_q_con,
2865       !    .             rain_con, snow_con, ibas_con, itop_con)
2866    ELSE IF (iflag_con.EQ.2) THEN
[3435]2867       CALL conflx(phys_tstep, paprs, pplay, t_seri, q_seri, &
[2469]2868            conv_t, conv_q, -evap, omega, &
2869            d_t_con, d_q_con, rain_con, snow_con, &
2870            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
2871            kcbot, kctop, kdtop, pmflxr, pmflxs)
2872       d_u_con = 0.
2873       d_v_con = 0.
[879]2874
[2469]2875       WHERE (rain_con < 0.) rain_con = 0.
2876       WHERE (snow_con < 0.) snow_con = 0.
2877       DO i = 1, klon
2878          ibas_con(i) = klev+1 - kcbot(i)
2879          itop_con(i) = klev+1 - kctop(i)
2880       ENDDO
2881    ELSE IF (iflag_con.GE.3) THEN
2882       ! nb of tracers for the KE convection:
2883       ! MAF la partie traceurs est faite dans phytrac
2884       ! on met ntra=1 pour limiter les appels mais on peut
2885       ! supprimer les calculs / ftra.
2886       ntra = 1
2887
2888       !=======================================================================
2889       !ajout pour la parametrisation des poches froides: calcul de
[2635]2890       !t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
[2692]2891       IF (iflag_wake>=1) THEN
2892         DO k=1,klev
2893            DO i=1,klon
2894                t_w(i,k) = t_seri(i,k) + (1-wake_s(i))*wake_deltat(i,k)
2895                q_w(i,k) = q_seri(i,k) + (1-wake_s(i))*wake_deltaq(i,k)
2896                t_x(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2897                q_x(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2898            ENDDO
2899         ENDDO
2900       ELSE
[4056]2901                t_w(:,:) = t_seri(:,:)
[2635]2902                q_w(:,:) = q_seri(:,:)
2903                t_x(:,:) = t_seri(:,:)
2904                q_x(:,:) = q_seri(:,:)
[2692]2905       ENDIF
[2469]2906       !
2907       !jyg<
2908       ! Perform dry adiabatic adjustment on wake profile
2909       ! The corresponding tendencies are added to the convective tendencies
2910       ! after the call to the convective scheme.
2911       IF (iflag_wake>=1) then
[2882]2912          IF (iflag_adjwk >= 1) THEN
[2469]2913             limbas(:) = 1
[2635]2914             CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
[2309]2915                  d_t_adjwk, d_q_adjwk)
[2638]2916             !
2917             DO k=1,klev
2918                DO i=1,klon
2919                   IF (wake_s(i) .GT. 1.e-3) THEN
2920                      t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
2921                      q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
2922                      d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
2923                      d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
2924                   ELSE
2925                      d_deltat_ajs_cv(i,k) = 0.
2926                      d_deltaq_ajs_cv(i,k) = 0.
2927                   ENDIF
2928                ENDDO
[2469]2929             ENDDO
[2882]2930             IF (iflag_adjwk == 2) THEN
2931               CALL add_wake_tend &
[3208]2932                 (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy)
[2882]2933             ENDIF  ! (iflag_adjwk == 2)
2934          ENDIF  ! (iflag_adjwk >= 1)
[2469]2935       ENDIF ! (iflag_wake>=1)
2936       !>jyg
2937       !
[2638]2938       
2939!!      print *,'physiq. q_w(1,k), q_x(1,k) ', &
2940!!             (k, q_w(1,k), q_x(1,k),k=1,25)
2941
[2513]2942!jyg<
[3435]2943       CALL alpale( debut, itap, phys_tstep, paprs, omega, t_seri,   &
[2513]2944                    alp_offset, it_wape_prescr,  wape_prescr, fip_prescr, &
2945                    ale_bl_prescr, alp_bl_prescr, &
2946                    wake_pe, wake_fip,  &
2947                    Ale_bl, Ale_bl_trig, Alp_bl, &
[2554]2948                    Ale, Alp , Ale_wake, Alp_wake)
[2513]2949!>jyg
2950!
[2469]2951       ! sb, oct02:
2952       ! Schema de convection modularise et vectorise:
2953       ! (driver commun aux versions 3 et 4)
2954       !
2955       IF (ok_cvl) THEN ! new driver for convectL
2956          !
2957          !jyg<
2958          ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2959          ! Calculate the upmost level of deep convection loops: k_upper_cv
2960          !  (near 22 km)
2961          k_upper_cv = klev
[3199]2962          !izero = klon/2+1/klon
2963          !DO k = klev,1,-1
2964          !   IF (pphi(izero,k) > 22.e4) k_upper_cv = k
2965          !ENDDO
2966          ! FH : nouveau calcul base sur un profil global sans quoi
2967          ! le modele etait sensible au decoupage de domaines
[2469]2968          DO k = klev,1,-1
[3199]2969             IF (-7*log(presnivs(k)/presnivs(1)) > 25.) k_upper_cv = k
[2469]2970          ENDDO
2971          IF (prt_level .ge. 5) THEN
2972             Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
2973                  k_upper_cv
2974          ENDIF
2975          !
2976          !>jyg
2977          IF (type_trac == 'repr') THEN
2978             nbtr_tmp=ntra
2979          ELSE
2980             nbtr_tmp=nbtr
[2692]2981          ENDIF
[2469]2982          !jyg   iflag_con est dans clesphys
2983          !c          CALL concvl (iflag_con,iflag_clos,
2984          CALL concvl (iflag_clos, &
[3435]2985               phys_tstep, paprs, pplay, k_upper_cv, t_x,q_x, &
[2635]2986               t_w,q_w,wake_s, &
[2469]2987               u_seri,v_seri,tr_seri,nbtr_tmp, &
2988               ALE,ALP, &
2989               sig1,w01, &
2990               d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2991               rain_con, snow_con, ibas_con, itop_con, sigd, &
[2824]2992               ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
[3496]2993               Ma,mipsh,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
[2469]2994               pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2995                                ! RomP >>>
2996                                !!     .        pmflxr,pmflxs,da,phi,mp,
2997                                !!     .        ftd,fqd,lalim_conv,wght_th)
[3496]2998               pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,qtaa,clw,elij, &
[2469]2999               ftd,fqd,lalim_conv,wght_th, &
3000               ev, ep,epmlmMm,eplaMm, &
[3496]3001               wdtrainA, wdtrainS, wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
[2481]3002               tau_cld_cv,coefw_cld_cv,epmax_diag)
[2630]3003
[2469]3004          ! RomP <<<
[619]3005
[2469]3006          !IM begin
3007          !       print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
3008          !    .dnwd0(1,1),ftd(1,1),fqd(1,1)
3009          !IM end
3010          !IM cf. FH
3011          clwcon0=qcondc
3012          pmfu(:,:)=upwd(:,:)+dnwd(:,:)
[3150]3013          !
3014          !jyg<
3015          ! If convective tendencies are too large, then call convection
3016          !  every time step
3017          cvpas = cvpas_0
3018          DO k=1,k_upper_cv
3019             DO i=1,klon
[3161]3020               IF (d_t_con(i,k) > 6.721 .AND. d_t_con(i,k) < 6.722 .AND.&
3021                   d_q_con(i,k) > -.0002171 .AND. d_q_con(i,k) < -.0002170) THEN
3022                     dtcon_multistep_max = 3.
3023                     dqcon_multistep_max = 0.02
3024               ENDIF
3025             ENDDO
3026          ENDDO
3027!
3028          DO k=1,k_upper_cv
3029             DO i=1,klon
[3150]3030!!               IF (abs(d_t_con(i,k)) > 0.24 .OR. &
3031!!                   abs(d_q_con(i,k)) > 2.e-2) THEN
3032               IF (abs(d_t_con(i,k)) > dtcon_multistep_max .OR. &
3033                   abs(d_q_con(i,k)) > dqcon_multistep_max) THEN
3034                 cvpas = 1
3035!!                 print *,'physiq1, i,k,d_t_con(i,k),d_q_con(i,k) ', &
3036!!                                   i,k,d_t_con(i,k),d_q_con(i,k)
3037               ENDIF
3038             ENDDO
3039          ENDDO
[3153]3040!!!   Ligne a ne surtout pas remettre sans avoir murement reflechi (jyg)
3041!!!          call bcast(cvpas)
3042!!!   ------------------------------------------------------------
[3150]3043          !>jyg
3044          !
[2692]3045          DO i = 1, klon
[3148]3046             IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+cvpas
[2692]3047          ENDDO
[2469]3048          !
3049          !jyg<
3050          !    Add the tendency due to the dry adjustment of the wake profile
3051          IF (iflag_wake>=1) THEN
[2882]3052            IF (iflag_adjwk == 2) THEN
3053              DO k=1,klev
3054                 DO i=1,klon
[3435]3055                    ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/phys_tstep
3056                    fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/phys_tstep
[2882]3057                    d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
3058                    d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
3059                 ENDDO
3060              ENDDO
3061            ENDIF  ! (iflag_adjwk = 2)
3062          ENDIF   ! (iflag_wake>=1)
[2469]3063          !>jyg
3064          !
3065       ELSE ! ok_cvl
[1412]3066
[2469]3067          ! MAF conema3 ne contient pas les traceurs
[3435]3068          CALL conema3 (phys_tstep, &
[2469]3069               paprs,pplay,t_seri,q_seri, &
3070               u_seri,v_seri,tr_seri,ntra, &
3071               sig1,w01, &
3072               d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
3073               rain_con, snow_con, ibas_con, itop_con, &
3074               upwd,dnwd,dnwd0,bas,top, &
3075               Ma,cape,tvp,rflag, &
3076               pbase &
3077               ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
3078               ,clwcon0)
[524]3079
[2469]3080       ENDIF ! ok_cvl
[524]3081
[2469]3082       !
3083       ! Correction precip
3084       rain_con = rain_con * cvl_corr
3085       snow_con = snow_con * cvl_corr
3086       !
[766]3087
[2469]3088       IF (.NOT. ok_gust) THEN
3089          do i = 1, klon
3090             wd(i)=0.0
3091          enddo
3092       ENDIF
[524]3093
[2469]3094       ! =================================================================== c
3095       ! Calcul des proprietes des nuages convectifs
3096       !
[524]3097
[2469]3098       !   calcul des proprietes des nuages convectifs
3099       clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
3100       IF (iflag_cld_cv == 0) THEN
[2692]3101          CALL clouds_gno &
[2469]3102               (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
3103       ELSE
[2692]3104          CALL clouds_bigauss &
[2469]3105               (klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
3106       ENDIF
[524]3107
[2205]3108
[2469]3109       ! =================================================================== c
[524]3110
[2469]3111       DO i = 1, klon
3112          itop_con(i) = min(max(itop_con(i),1),klev)
3113          ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
3114       ENDDO
[1428]3115
[2469]3116       DO i = 1, klon
[4056]3117          ! C Risi modif: pour éviter pb de dépassement d'indice dans les cas
3118          ! où i n'est pas un point convectif et donc ibas_con(i)=0
3119          ! c'est un pb indépendant des isotopes
3120          if (ibas_con(i) > 0) then
3121             ema_pcb(i)  = paprs(i,ibas_con(i))
3122          else
3123             ema_pcb(i)  = 0.0
3124          endif
[2469]3125       ENDDO
3126       DO i = 1, klon
3127          ! L'idicage de itop_con peut cacher un pb potentiel
3128          ! FH sous la dictee de JYG, CR
3129          ema_pct(i)  = paprs(i,itop_con(i)+1)
[879]3130
[2692]3131          IF (itop_con(i).gt.klev-3) THEN
3132             IF (prt_level >= 9) THEN
[2469]3133                write(lunout,*)'La convection monte trop haut '
3134                write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
[2692]3135             ENDIF
3136          ENDIF
[2469]3137       ENDDO
3138    ELSE IF (iflag_con.eq.0) THEN
3139       write(lunout,*) 'On n appelle pas la convection'
3140       clwcon0=0.
3141       rnebcon0=0.
3142       d_t_con=0.
3143       d_q_con=0.
3144       d_u_con=0.
3145       d_v_con=0.
3146       rain_con=0.
3147       snow_con=0.
3148       bas=1
3149       top=1
3150    ELSE
3151       WRITE(lunout,*) "iflag_con non-prevu", iflag_con
[2692]3152       CALL abort_physic("physiq", "", 1)
[2469]3153    ENDIF
[524]3154
[2469]3155    !     CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
3156    !    .              d_u_con, d_v_con)
[524]3157
[2730]3158!jyg    Reinitialize proba_notrig and itapcv when convection has been called
3159    proba_notrig(:) = 1.
[2707]3160    itapcv = 0
[3150]3161    ENDIF !  (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itapcv,cvpas_0).EQ.0)
[2730]3162!
[2707]3163    itapcv = itapcv+1
[3153]3164    !
3165    ! Compter les steps ou cvpas=1
3166    IF (cvpas == 1) THEN
3167      Ncvpaseq1 = Ncvpaseq1+1
3168    ENDIF
3169    IF (mod(itap,1000) == 0) THEN
3170      print *,' physiq, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
3171    ENDIF
[2707]3172
[2812]3173!!!jyg  Appel diagnostique a add_phys_tend pour tester la conservation de
3174!!!     l'energie dans les courants satures.
3175!!    d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
3176!!    d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
3177!!    dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
3178!!    CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat,   &
3179!!                     dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,& 
3180!!                     itap, 1)
3181!!    call prt_enerbil('convection_sat',itap)
3182!!
3183!!
[2469]3184    CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
[2812]3185         'convection',abortphy,flag_inhib_tend,itap,0)
[3461]3186    CALL prt_enerbil('convection',itap)
[2235]3187
[2469]3188    !-------------------------------------------------------------------------
[766]3189
[2692]3190    IF (mydebug) THEN
3191       CALL writefield_phy('u_seri',u_seri,nbp_lev)
3192       CALL writefield_phy('v_seri',v_seri,nbp_lev)
3193       CALL writefield_phy('t_seri',t_seri,nbp_lev)
3194       CALL writefield_phy('q_seri',q_seri,nbp_lev)
3195    ENDIF
[766]3196
[2469]3197    IF (check) THEN
3198       za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3199       WRITE(lunout,*)"aprescon=", za
3200       zx_t = 0.0
3201       za = 0.0
3202       DO i = 1, klon
3203          za = za + cell_area(i)/REAL(klon)
3204          zx_t = zx_t + (rain_con(i)+ &
3205               snow_con(i))*cell_area(i)/REAL(klon)
3206       ENDDO
[3435]3207       zx_t = zx_t/za*phys_tstep
[2469]3208       WRITE(lunout,*)"Precip=", zx_t
3209    ENDIF
3210    IF (zx_ajustq) THEN
3211       DO i = 1, klon
3212          z_apres(i) = 0.0
3213       ENDDO
3214       DO k = 1, klev
3215          DO i = 1, klon
3216             z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
3217                  *(paprs(i,k)-paprs(i,k+1))/RG
3218          ENDDO
3219       ENDDO
3220       DO i = 1, klon
[3435]3221          z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*phys_tstep) &
[2469]3222               /z_apres(i)
3223       ENDDO
3224       DO k = 1, klev
3225          DO i = 1, klon
3226             IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
3227                  z_factor(i).LT.(1.0-1.0E-08)) THEN
3228                q_seri(i,k) = q_seri(i,k) * z_factor(i)
3229             ENDIF
3230          ENDDO
3231       ENDDO
3232    ENDIF
3233    zx_ajustq=.FALSE.
[879]3234
[2469]3235    !
3236    !==========================================================================
3237    !RR:Evolution de la poche froide: on ne fait pas de separation wake/env
3238    !pour la couche limite diffuse pour l instant
3239    !
3240    !
3241    ! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
3242    ! il faut rajouter cette tendance calcul\'ee hors des poches
3243    ! froides
3244    !
[2692]3245    IF (iflag_wake>=1) THEN
[2707]3246       !
3247       !
[2730]3248       ! Call wakes every "wkpas" step
3249       !
3250       IF (MOD(itapwk,wkpas).EQ.0) THEN
3251          !
3252          DO k=1,klev
[2469]3253             DO i=1,klon
[2730]3254                dt_dwn(i,k)  = ftd(i,k)
3255                dq_dwn(i,k)  = fqd(i,k)
3256                M_dwn(i,k)   = dnwd0(i,k)
3257                M_up(i,k)    = upwd(i,k)
[3435]3258                dt_a(i,k)    = d_t_con(i,k)/phys_tstep - ftd(i,k)
3259                dq_a(i,k)    = d_q_con(i,k)/phys_tstep - fqd(i,k)
[2469]3260             ENDDO
3261          ENDDO
[2730]3262         
3263          IF (iflag_wake==2) THEN
3264             ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3265             DO k = 1,klev
3266                dt_dwn(:,k)= dt_dwn(:,k)+ &
[3435]3267                     ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/phys_tstep
[2730]3268                dq_dwn(:,k)= dq_dwn(:,k)+ &
[3435]3269                     ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/phys_tstep
[2730]3270             ENDDO
3271          ELSEIF (iflag_wake==3) THEN
3272             ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3273             DO k = 1,klev
3274                DO i=1,klon
3275                   IF (rneb(i,k)==0.) THEN
3276                      ! On ne tient compte des tendances qu'en dehors des
3277                      ! nuages (c'est-\`a-dire a priri dans une region ou
3278                      ! l'eau se reevapore).
3279                      dt_dwn(i,k)= dt_dwn(i,k)+ &
[3435]3280                           ok_wk_lsp(i)*d_t_lsc(i,k)/phys_tstep
[2730]3281                      dq_dwn(i,k)= dq_dwn(i,k)+ &
[3435]3282                           ok_wk_lsp(i)*d_q_lsc(i,k)/phys_tstep
[2730]3283                   ENDIF
3284                ENDDO
3285             ENDDO
3286          ENDIF
3287         
3288          !
3289          !calcul caracteristiques de la poche froide
[3435]3290          CALL calWAKE (iflag_wake_tend, paprs, pplay, phys_tstep, &
[2730]3291               t_seri, q_seri, omega,  &
3292               dt_dwn, dq_dwn, M_dwn, M_up,  &
[3208]3293               dt_a, dq_a, cv_gen,  &
3294               sigd, cin,  &
3295               wake_deltat, wake_deltaq, wake_s, awake_dens, wake_dens,  &
[2730]3296               wake_dth, wake_h,  &
[3000]3297!!               wake_pe, wake_fip, wake_gfl,  &
3298               wake_pe, wake_fip_0, wake_gfl,  &   !! jyg
[2730]3299               d_t_wake, d_q_wake,  &
3300               wake_k, t_x, q_x,  &
3301               wake_omgbdth, wake_dp_omgb,  &
3302               wake_dtKE, wake_dqKE,  &
3303               wake_omg, wake_dp_deltomg,  &
3304               wake_spread, wake_Cstar, d_deltat_wk_gw,  &
[3208]3305               d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk)
[2730]3306          !
3307          !jyg    Reinitialize itapwk when wakes have been called
3308          itapwk = 0
3309       ENDIF !  (MOD(itapwk,wkpas).EQ.0)
[2469]3310       !
[2730]3311       itapwk = itapwk+1
[2469]3312       !
3313       !-----------------------------------------------------------------------
3314       ! ajout des tendances des poches froides
3315       CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
[2812]3316            abortphy,flag_inhib_tend,itap,0)
[3461]3317       CALL prt_enerbil('wake',itap)
[2469]3318       !------------------------------------------------------------------------
[879]3319
[2730]3320       ! Increment Wake state variables
[2635]3321       IF (iflag_wake_tend .GT. 0.) THEN
3322
3323         CALL add_wake_tend &
[3208]3324            (d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk, wake_k, &
[2635]3325             'wake', abortphy)
[3461]3326          CALL prt_enerbil('wake',itap)
[2635]3327       ENDIF   ! (iflag_wake_tend .GT. 0.)
[3179]3328       !
3329       IF (prt_level .GE. 10) THEN
3330         print *,' physiq, after calwake, wake_s: ',wake_s(:)
3331         print *,' physiq, after calwake, wake_deltat: ',wake_deltat(:,1)
3332         print *,' physiq, after calwake, wake_deltaq: ',wake_deltaq(:,1)
3333       ENDIF
[2635]3334
[3000]3335       IF (iflag_alp_wk_cond .GT. 0.) THEN
3336
[3435]3337         CALL alpale_wk(phys_tstep, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
[3000]3338                        wake_fip)
3339       ELSE
3340         wake_fip(:) = wake_fip_0(:)
3341       ENDIF   ! (iflag_alp_wk_cond .GT. 0.)
3342
[2692]3343    ENDIF  ! (iflag_wake>=1)
[2469]3344    !
3345    !===================================================================
3346    ! Convection seche (thermiques ou ajustement)
3347    !===================================================================
3348    !
[2692]3349    CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
[2469]3350         ,seuil_inversion,weak_inversion,dthmin)
[878]3351
3352
3353
[2469]3354    d_t_ajsb(:,:)=0.
3355    d_q_ajsb(:,:)=0.
3356    d_t_ajs(:,:)=0.
3357    d_u_ajs(:,:)=0.
3358    d_v_ajs(:,:)=0.
3359    d_q_ajs(:,:)=0.
3360    clwcon0th(:,:)=0.
3361    !
3362    !      fm_therm(:,:)=0.
3363    !      entr_therm(:,:)=0.
3364    !      detr_therm(:,:)=0.
3365    !
[2692]3366    IF (prt_level>9) WRITE(lunout,*) &
[2469]3367         'AVANT LA CONVECTION SECHE , iflag_thermals=' &
3368         ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
[2692]3369    IF (iflag_thermals<0) THEN
[2469]3370       !  Rien
3371       !  ====
[2692]3372       IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
[541]3373
[878]3374
[2692]3375    ELSE
[878]3376
[2469]3377       !  Thermiques
3378       !  ==========
[2692]3379       IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
[2469]3380            ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
[878]3381
3382
[2469]3383       !cc nrlmd le 10/04/2012
3384       DO k=1,klev+1
3385          DO i=1,klon
3386             pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
3387             pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
3388             pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
3389             pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
[2159]3390          ENDDO
[2469]3391       ENDDO
3392       !cc fin nrlmd le 10/04/2012
[1403]3393
[2692]3394       IF (iflag_thermals>=1) THEN
[2469]3395          !jyg<
[2852]3396!!       IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3397       IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
[2469]3398             !  Appel des thermiques avec les profils exterieurs aux poches
3399             DO k=1,klev
3400                DO i=1,klon
3401                   t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
3402                   q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
[2606]3403                   u_therm(i,k) = u_seri(i,k)
3404                   v_therm(i,k) = v_seri(i,k)
[2469]3405                ENDDO
3406             ENDDO
3407          ELSE
3408             !  Appel des thermiques avec les profils moyens
3409             DO k=1,klev
3410                DO i=1,klon
3411                   t_therm(i,k) = t_seri(i,k)
3412                   q_therm(i,k) = q_seri(i,k)
[2606]3413                   u_therm(i,k) = u_seri(i,k)
3414                   v_therm(i,k) = v_seri(i,k)
[2469]3415                ENDDO
3416             ENDDO
3417          ENDIF
3418          !>jyg
[2692]3419          CALL calltherm(pdtphys &
[2469]3420               ,pplay,paprs,pphi,weak_inversion &
[2606]3421                        ! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
3422               ,u_therm,v_therm,t_therm,q_therm,zqsat,debut &  !jyg
[2469]3423               ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
3424               ,fm_therm,entr_therm,detr_therm &
3425               ,zqasc,clwcon0th,lmax_th,ratqscth &
3426               ,ratqsdiff,zqsatth &
3427                                !on rajoute ale et alp, et les
3428                                !caracteristiques de la couche alim
3429               ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
3430               ,ztv,zpspsk,ztla,zthl &
3431                                !cc nrlmd le 10/04/2012
3432               ,pbl_tke_input,pctsrf,omega,cell_area &
3433               ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
3434               ,n2,s2,ale_bl_stat &
3435               ,therm_tke_max,env_tke_max &
3436               ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
3437               ,alp_bl_conv,alp_bl_stat &
3438                                !cc fin nrlmd le 10/04/2012
3439               ,zqla,ztva )
3440          !
3441          !jyg<
[2852]3442!!jyg          IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3443          IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
[2469]3444             !  Si les thermiques ne sont presents que hors des
3445             !  poches, la tendance moyenne associ\'ee doit etre
3446             !  multipliee par la fraction surfacique qu'ils couvrent.
3447             DO k=1,klev
3448                DO i=1,klon
3449                   !
[2635]3450                   d_deltat_the(i,k) = - d_t_ajs(i,k)
3451                   d_deltaq_the(i,k) = - d_q_ajs(i,k)
[2469]3452                   !
3453                   d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
3454                   d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
3455                   d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
3456                   d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
3457                   !
3458                ENDDO
3459             ENDDO
[2606]3460          !
[3180]3461             IF (ok_bug_split_th) THEN
3462               CALL add_wake_tend &
[3208]3463                   (d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy)
[3180]3464             ELSE
3465               CALL add_wake_tend &
[3208]3466                   (d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wake_k, 'the', abortphy)
[3180]3467             ENDIF
[3461]3468             CALL prt_enerbil('the',itap)
[2638]3469          !
[2852]3470          ENDIF  ! (mod(iflag_pbl_split/10,10) .GE. 1)
[2638]3471          !
[2606]3472          CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,  &
[2812]3473                             dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
[3461]3474          CALL prt_enerbil('thermals',itap)
[2606]3475          !
[2513]3476!
[3435]3477          CALL alpale_th( phys_tstep, lmax_th, t_seri, cell_area,  &
[2513]3478                          cin, s2, n2,  &
3479                          ale_bl_trig, ale_bl_stat, ale_bl,  &
[2556]3480                          alp_bl, alp_bl_stat, &
[3208]3481                          proba_notrig, random_notrig, cv_gen)
[2635]3482          !>jyg
[1638]3483
[2554]3484          ! ------------------------------------------------------------------
3485          ! Transport de la TKE par les panaches thermiques.
3486          ! FH : 2010/02/01
3487          !     if (iflag_pbl.eq.10) then
3488          !     call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
3489          !    s           rg,paprs,pbl_tke)
3490          !     endif
3491          ! -------------------------------------------------------------------
3492
[2692]3493          DO i=1,klon
[2469]3494             !           zmax_th(i)=pphi(i,lmax_th(i))/rg
3495             !CR:04/05/12:correction calcul zmax
3496             zmax_th(i)=zmax0(i)
[2692]3497          ENDDO
[1507]3498
[2692]3499       ENDIF
[878]3500
[2469]3501       !  Ajustement sec
3502       !  ==============
[878]3503
[2469]3504       ! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3505       ! a partir du sommet des thermiques.
3506       ! Dans le cas contraire, on demarre au niveau 1.
[878]3507
[2692]3508       IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
[878]3509
[2692]3510          IF (iflag_thermals.eq.0) THEN
3511             IF (prt_level>9) WRITE(lunout,*)'ajsec'
[2469]3512             limbas(:)=1
[2692]3513          ELSE
[2469]3514             limbas(:)=lmax_th(:)
[2692]3515          ENDIF
[878]3516
[2469]3517          ! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3518          ! pour des test de convergence numerique.
3519          ! Le nouveau ajsec est a priori mieux, meme pour le cas
3520          ! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3521          ! non nulles numeriquement pour des mailles non concernees.
[878]3522
[2692]3523          IF (iflag_thermals==0) THEN
[2469]3524             ! Calling adjustment alone (but not the thermal plume model)
3525             CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3526                  , d_t_ajsb, d_q_ajsb)
[2692]3527          ELSE IF (iflag_thermals>0) THEN
[2469]3528             ! Calling adjustment above the top of thermal plumes
3529             CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3530                  , d_t_ajsb, d_q_ajsb)
[2692]3531          ENDIF
[878]3532
[2469]3533          !--------------------------------------------------------------------
3534          ! ajout des tendances de l'ajustement sec ou des thermiques
3535          CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
[2812]3536               'ajsb',abortphy,flag_inhib_tend,itap,0)
[3461]3537          CALL prt_enerbil('ajsb',itap)
[2469]3538          d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3539          d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
[904]3540
[2469]3541          !---------------------------------------------------------------------
[878]3542
[2692]3543       ENDIF
[524]3544
[2692]3545    ENDIF
[2469]3546    !
3547    !===================================================================
3548    ! Computation of ratqs, the width (normalized) of the subrid scale
3549    ! water distribution
[4009]3550
3551    tke_dissip_ave(:,:)=0.
3552    l_mix_ave(:,:)=0.
3553    wprime_ave(:,:)=0.
3554
3555    DO nsrf = 1, nbsrf
3556       DO i = 1, klon
3557          tke_dissip_ave(i,:) = tke_dissip_ave(i,:) + tke_dissip(i,:,nsrf)*pctsrf(i,nsrf)
3558          l_mix_ave(i,:) = l_mix_ave(i,:) + l_mix(i,:,nsrf)*pctsrf(i,nsrf)
3559          wprime_ave(i,:) = wprime_ave(i,:) + wprime(i,:,nsrf)*pctsrf(i,nsrf)
3560       ENDDO
3561    ENDDO
3562
[2469]3563    CALL  calcratqs(klon,klev,prt_level,lunout,        &
3564         iflag_ratqs,iflag_con,iflag_cld_th,pdtphys,  &
[2534]3565         ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
[3856]3566         tau_ratqs,fact_cldcon,wake_s, wake_deltaq,   &
[2469]3567         ptconv,ptconvth,clwcon0th, rnebcon0th,     &
[4009]3568         paprs,pplay,t_seri,q_seri, qtc_cv, sigt_cv, zqsat, &
3569         pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave,t2m,q2m,fm_therm, &
[3856]3570         ratqs,ratqsc,ratqs_inter)
[1032]3571
[2469]3572    !
3573    ! Appeler le processus de condensation a grande echelle
3574    ! et le processus de precipitation
3575    !-------------------------------------------------------------------------
3576    IF (prt_level .GE.10) THEN
3577       print *,'itap, ->fisrtilp ',itap
3578    ENDIF
[4056]3579    !
[3999]3580
3581    picefra(:,:)=0.
3582
3583    IF (ok_new_lscp) THEN
3584
[4062]3585    !--mise à jour de flight_m et flight_h2o dans leur module
3586    IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
3587      CALL airplane(debut,pphis,pplay,paprs,t_seri)
3588    ENDIF
[4059]3589
3590    CALL lscp(phys_tstep,missing_val,paprs,pplay, &
[3999]3591         t_seri, q_seri,ptconv,ratqs, &
[4059]3592         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneb_seri, &
3593         cldliq, picefra, rain_lsc, snow_lsc, &
[3999]3594         pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3595         frac_impa, frac_nucl, beta_prec_fisrt, &
3596         prfl, psfl, rhcl,  &
3597         zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
[4062]3598         iflag_ice_thermo, ok_ice_sursat)
[3999]3599
3600    ELSE
[4059]3601
[3435]3602    CALL fisrtilp(phys_tstep,paprs,pplay, &
[2469]3603         t_seri, q_seri,ptconv,ratqs, &
3604         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
3605         rain_lsc, snow_lsc, &
3606         pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3607         frac_impa, frac_nucl, beta_prec_fisrt, &
3608         prfl, psfl, rhcl,  &
3609         zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3610         iflag_ice_thermo)
[4059]3611
[3999]3612    ENDIF
[4056]3613    !
[2469]3614    WHERE (rain_lsc < 0) rain_lsc = 0.
3615    WHERE (snow_lsc < 0) snow_lsc = 0.
[766]3616
[2799]3617!+JLD
3618!    write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
3619!        & ,((rcw-rcpd)*rain_lsc + (rcs-rcpd)*snow_lsc)*t_seri(1,1)-rlvtt*rain_lsc+rlstt*snow_lsc &
3620!        & ,rain_lsc,snow_lsc
3621!    write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
3622!-JLD
[2469]3623    CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
[2812]3624         'lsc',abortphy,flag_inhib_tend,itap,0)
[3461]3625    CALL prt_enerbil('lsc',itap)
[2613]3626    rain_num(:)=0.
[2657]3627    DO k = 1, klev
[2613]3628       DO i = 1, klon
3629          IF (ql_seri(i,k)>oliqmax) THEN
3630             rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
3631             ql_seri(i,k)=oliqmax
3632          ENDIF
3633       ENDDO
3634    ENDDO
[4098]3635    IF (nqo >= 3) THEN
[2657]3636    DO k = 1, klev
3637       DO i = 1, klon
3638          IF (qs_seri(i,k)>oicemax) THEN
3639             rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
3640             qs_seri(i,k)=oicemax
3641          ENDIF
3642       ENDDO
3643    ENDDO
3644    ENDIF
[2613]3645
[2524]3646    !---------------------------------------------------------------------------
[2469]3647    DO k = 1, klev
3648       DO i = 1, klon
3649          cldfra(i,k) = rneb(i,k)
3650          !CR: a quoi ca sert? Faut-il ajouter qs_seri?
3651          IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3652       ENDDO
3653    ENDDO
3654    IF (check) THEN
3655       za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3656       WRITE(lunout,*)"apresilp=", za
3657       zx_t = 0.0
3658       za = 0.0
3659       DO i = 1, klon
3660          za = za + cell_area(i)/REAL(klon)
3661          zx_t = zx_t + (rain_lsc(i) &
3662               + snow_lsc(i))*cell_area(i)/REAL(klon)
3663       ENDDO
[3435]3664       zx_t = zx_t/za*phys_tstep
[2469]3665       WRITE(lunout,*)"Precip=", zx_t
3666    ENDIF
[766]3667
[2692]3668    IF (mydebug) THEN
3669       CALL writefield_phy('u_seri',u_seri,nbp_lev)
3670       CALL writefield_phy('v_seri',v_seri,nbp_lev)
3671       CALL writefield_phy('t_seri',t_seri,nbp_lev)
3672       CALL writefield_phy('q_seri',q_seri,nbp_lev)
3673    ENDIF
[524]3674
[2469]3675    !
3676    !-------------------------------------------------------------------
3677    !  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3678    !-------------------------------------------------------------------
[524]3679
[2469]3680    ! 1. NUAGES CONVECTIFS
3681    !
3682    !IM cf FH
3683    !     IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
3684    IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
3685       snow_tiedtke=0.
3686       !     print*,'avant calcul de la pseudo precip '
3687       !     print*,'iflag_cld_th',iflag_cld_th
[2692]3688       IF (iflag_cld_th.eq.-1) THEN
[2469]3689          rain_tiedtke=rain_con
[2692]3690       ELSE
[2469]3691          !       print*,'calcul de la pseudo precip '
3692          rain_tiedtke=0.
3693          !         print*,'calcul de la pseudo precip 0'
[2692]3694          DO k=1,klev
3695             DO i=1,klon
3696                IF (d_q_con(i,k).lt.0.) THEN
[2469]3697                   rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3698                        *(paprs(i,k)-paprs(i,k+1))/rg
[2692]3699                ENDIF
3700             ENDDO
3701          ENDDO
3702       ENDIF
[2469]3703       !
3704       !     call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3705       !
[524]3706
[2469]3707       ! Nuages diagnostiques pour Tiedtke
3708       CALL diagcld1(paprs,pplay, &
3709                                !IM cf FH. rain_con,snow_con,ibas_con,itop_con,
3710            rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3711            diafra,dialiq)
3712       DO k = 1, klev
3713          DO i = 1, klon
3714             IF (diafra(i,k).GT.cldfra(i,k)) THEN
3715                cldliq(i,k) = dialiq(i,k)
3716                cldfra(i,k) = diafra(i,k)
3717             ENDIF
3718          ENDDO
3719       ENDDO
[524]3720
[2469]3721    ELSE IF (iflag_cld_th.ge.3) THEN
3722       !  On prend pour les nuages convectifs le max du calcul de la
3723       !  convection et du calcul du pas de temps precedent diminue d'un facteur
3724       !  facttemps
3725       facteur = pdtphys *facttemps
[2692]3726       DO k=1,klev
3727          DO i=1,klon
[2469]3728             rnebcon(i,k)=rnebcon(i,k)*facteur
[2692]3729             IF (rnebcon0(i,k)*clwcon0(i,k).GT.rnebcon(i,k)*clwcon(i,k)) THEN
[2469]3730                rnebcon(i,k)=rnebcon0(i,k)
3731                clwcon(i,k)=clwcon0(i,k)
[2692]3732             ENDIF
3733          ENDDO
3734       ENDDO
[2469]3735
3736       !   On prend la somme des fractions nuageuses et des contenus en eau
[524]3737
[2692]3738       IF (iflag_cld_th>=5) THEN
[1411]3739
[2692]3740          DO k=1,klev
[2469]3741             ptconvth(:,k)=fm_therm(:,k+1)>0.
[2692]3742          ENDDO
[1496]3743
[2692]3744          IF (iflag_coupl==4) THEN
[1496]3745
[2469]3746             ! Dans le cas iflag_coupl==4, on prend la somme des convertures
3747             ! convectives et lsc dans la partie des thermiques
3748             ! Le controle par iflag_coupl est peut etre provisoire.
[2692]3749             DO k=1,klev
3750                DO i=1,klon
3751                   IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
[2469]3752                      cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3753                      cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
[2692]3754                   ELSE IF (ptconv(i,k)) THEN
[2469]3755                      cldfra(i,k)=rnebcon(i,k)
3756                      cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
[2692]3757                   ENDIF
3758                ENDDO
3759             ENDDO
[1496]3760
[2692]3761          ELSE IF (iflag_coupl==5) THEN
3762             DO k=1,klev
3763                DO i=1,klon
[2469]3764                   cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3765                   cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
[2692]3766                ENDDO
3767             ENDDO
[1525]3768
[2692]3769          ELSE
[1525]3770
[2469]3771             ! Si on est sur un point touche par la convection
3772             ! profonde et pas par les thermiques, on prend la
3773             ! couverture nuageuse et l'eau nuageuse de la convection
3774             ! profonde.
[1411]3775
[2469]3776             !IM/FH: 2011/02/23
3777             ! definition des points sur lesquels ls thermiques sont actifs
[1496]3778
[2692]3779             DO k=1,klev
3780                DO i=1,klon
3781                   IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
[2469]3782                      cldfra(i,k)=rnebcon(i,k)
3783                      cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
[2692]3784                   ENDIF
3785                ENDDO
3786             ENDDO
[1496]3787
[2692]3788          ENDIF
[1496]3789
[2692]3790       ELSE
[1496]3791
[2469]3792          ! Ancienne version
3793          cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3794          cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
[2692]3795       ENDIF
[1411]3796
[2469]3797    ENDIF
[1507]3798
[2469]3799    !     plulsc(:)=0.
3800    !     do k=1,klev,-1
3801    !        do i=1,klon
3802    !              zzz=prfl(:,k)+psfl(:,k)
3803    !           if (.not.ptconvth.zzz.gt.0.)
3804    !        enddo prfl, psfl,
3805    !     enddo
3806    !
3807    ! 2. NUAGES STARTIFORMES
3808    !
3809    IF (ok_stratus) THEN
3810       CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3811       DO k = 1, klev
3812          DO i = 1, klon
3813             IF (diafra(i,k).GT.cldfra(i,k)) THEN
3814                cldliq(i,k) = dialiq(i,k)
3815                cldfra(i,k) = diafra(i,k)
3816             ENDIF
3817          ENDDO
3818       ENDDO
3819    ENDIF
3820    !
3821    ! Precipitation totale
3822    !
3823    DO i = 1, klon
3824       rain_fall(i) = rain_con(i) + rain_lsc(i)
3825       snow_fall(i) = snow_con(i) + snow_lsc(i)
3826    ENDDO
3827    !
3828    ! Calculer l'humidite relative pour diagnostique
3829    !
3830    DO k = 1, klev
3831       DO i = 1, klon
3832          zx_t = t_seri(i,k)
3833          IF (thermcep) THEN
3834             !!           if (iflag_ice_thermo.eq.0) then                 !jyg
3835             zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3836             !!           else                                            !jyg
3837             !!           zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t))      !jyg
3838             !!           endif                                           !jyg
3839             zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3840             zx_qs  = MIN(0.5,zx_qs)
3841             zcor   = 1./(1.-retv*zx_qs)
3842             zx_qs  = zx_qs*zcor
3843          ELSE
3844             !!           IF (zx_t.LT.t_coup) THEN             !jyg
3845             IF (zx_t.LT.rtt) THEN                  !jyg
3846                zx_qs = qsats(zx_t)/pplay(i,k)
3847             ELSE
3848                zx_qs = qsatl(zx_t)/pplay(i,k)
3849             ENDIF
3850          ENDIF
3851          zx_rh(i,k) = q_seri(i,k)/zx_qs
[3784]3852            IF (iflag_ice_thermo .GT. 0) THEN
[3780]3853          zx_rhl(i,k) = q_seri(i,k)/(qsatl(zx_t)/pplay(i,k))
3854          zx_rhi(i,k) = q_seri(i,k)/(qsats(zx_t)/pplay(i,k))
[3784]3855            ENDIF
[2469]3856          zqsat(i,k)=zx_qs
3857       ENDDO
3858    ENDDO
[782]3859
[2469]3860    !IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3861    !   equivalente a 2m (tpote) pour diagnostique
3862    !
3863    DO i = 1, klon
3864       tpot(i)=zt2m(i)*(100000./paprs(i,1))**RKAPPA
3865       IF (thermcep) THEN
3866          IF(zt2m(i).LT.RTT) then
3867             Lheat=RLSTT
3868          ELSE
3869             Lheat=RLVTT
3870          ENDIF
3871       ELSE
3872          IF (zt2m(i).LT.RTT) THEN
3873             Lheat=RLSTT
3874          ELSE
3875             Lheat=RLVTT
3876          ENDIF
3877       ENDIF
3878       tpote(i) = tpot(i)*      &
3879            EXP((Lheat *qsat2m(i))/(RCPD*zt2m(i)))
3880    ENDDO
[524]3881
[3870]3882    IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN      ! ModThL
[524]3883#ifdef INCA
[2469]3884       CALL VTe(VTphysiq)
3885       CALL VTb(VTinca)
3886       calday = REAL(days_elapsed + 1) + jH_cur
[524]3887
[3435]3888       CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
[3563]3889       CALL AEROSOL_METEO_CALC( &
3890            calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3891            prfl,psfl,pctsrf,cell_area, &
3892            latitude_deg,longitude_deg,u10m,v10m)
[524]3893
[2469]3894       zxsnow_dummy(:) = 0.0
[625]3895
[2469]3896       CALL chemhook_begin (calday, &
3897            days_elapsed+1, &
3898            jH_cur, &
3899            pctsrf(1,1), &
3900            latitude_deg, &
3901            longitude_deg, &
3902            cell_area, &
3903            paprs, &
3904            pplay, &
3905            coefh(1:klon,1:klev,is_ave), &
3906            pphi, &
3907            t_seri, &
3908            u, &
3909            v, &
[3613]3910            rot, &
[2469]3911            wo(:, :, 1), &
3912            q_seri, &
3913            zxtsol, &
[3613]3914            zt2m, &
[2469]3915            zxsnow_dummy, &
3916            solsw, &
3917            albsol1, &
3918            rain_fall, &
3919            snow_fall, &
3920            itop_con, &
3921            ibas_con, &
3922            cldfra, &
3923            nbp_lon, &
3924            nbp_lat-1, &
[3872]3925            tr_seri(:,:,1+nqCO2:nbtr), &
[2469]3926            ftsol, &
3927            paprs, &
3928            cdragh, &
3929            cdragm, &
3930            pctsrf, &
3931            pdtphys, &
3932            itap)
[616]3933
[2469]3934       CALL VTe(VTinca)
3935       CALL VTb(VTphysiq)
[3865]3936#endif
3937    ENDIF !type_trac = inca or inco
[3666]3938    IF (type_trac == 'repr') THEN
3939#ifdef REPROBUS
3940    !CALL chemtime_rep(itap+itau_phy-1, date0, dtime, itap)
3941    CALL chemtime_rep(itap+itau_phy-1, date0, phys_tstep, itap)
3942#endif
3943    ENDIF
[2618]3944
[2469]3945    !
[2618]3946    ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3947    !
3948    IF (MOD(itaprad,radpas).EQ.0) THEN
[959]3949
[2618]3950       !
3951       !jq - introduce the aerosol direct and first indirect radiative forcings
3952       !jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
[2738]3953       IF (flag_aerosol .GT. 0) THEN
[2618]3954          IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3955             IF (.NOT. aerosol_couple) THEN
3956                !
3957                CALL readaerosol_optic( &
[3630]3958                     debut, flag_aerosol, itap, jD_cur-jD_ref, &
[2618]3959                     pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
3960                     mass_solu_aero, mass_solu_aero_pi,  &
3961                     tau_aero, piz_aero, cg_aero,  &
3962                     tausum_aero, tau3d_aero)
3963             ENDIF
3964          ELSE                       ! RRTM radiation
3965             IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
3966                abort_message='config_inca=aero et rrtm=1 impossible'
[2692]3967                CALL abort_physic(modname,abort_message,1)
[2618]3968             ELSE
3969                !
3970#ifdef CPP_RRTM
3971                IF (NSW.EQ.6) THEN
[2738]3972                   !--new aerosol properties SW and LW
[2618]3973                   !
[2753]3974#ifdef CPP_Dust
3975                   !--SPL aerosol model
3976                   CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
3977                        tr_seri, mass_solu_aero, mass_solu_aero_pi,  &
3978                        tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,  &
3979                        tausum_aero, tau3d_aero)
3980#else
3981                   !--climatologies or INCA aerosols
[3480]3982                   CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
[3630]3983                        flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
[2618]3984                        pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
3985                        tr_seri, mass_solu_aero, mass_solu_aero_pi,  &
3986                        tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,  &
[2854]3987                        tausum_aero, drytausum_aero, tau3d_aero)
[2753]3988#endif
[3274]3989
3990                   IF (flag_aerosol .EQ. 7) THEN
3991                      CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg,  &
3992                                   tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
3993                   ENDIF
3994
[2738]3995                   !
[2618]3996                ELSE IF (NSW.EQ.2) THEN
3997                   !--for now we use the old aerosol properties
3998                   !
3999                   CALL readaerosol_optic( &
[3630]4000                        debut, flag_aerosol, itap, jD_cur-jD_ref, &
[2618]4001                        pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
4002                        mass_solu_aero, mass_solu_aero_pi,  &
4003                        tau_aero, piz_aero, cg_aero,  &
4004                        tausum_aero, tau3d_aero)
4005                   !
4006                   !--natural aerosols
4007                   tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
4008                   piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
4009                   cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
4010                   !--all aerosols
4011                   tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
4012                   piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
4013                   cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
[2738]4014                   !
4015                   !--no LW optics
4016                   tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4017                   !
[2618]4018                ELSE
4019                   abort_message='Only NSW=2 or 6 are possible with ' &
4020                        // 'aerosols and iflag_rrtm=1'
[2692]4021                   CALL abort_physic(modname,abort_message,1)
[2618]4022                ENDIF
4023#else
4024                abort_message='You should compile with -rrtm if running ' &
4025                     // 'with iflag_rrtm=1'
[2692]4026                CALL abort_physic(modname,abort_message,1)
[2618]4027#endif
4028                !
4029             ENDIF
4030          ENDIF
[2738]4031       ELSE   !--flag_aerosol = 0
[2618]4032          tausum_aero(:,:,:) = 0.
[2854]4033          drytausum_aero(:,:) = 0.
[2640]4034          mass_solu_aero(:,:) = 0.
4035          mass_solu_aero_pi(:,:) = 0.
[2618]4036          IF (iflag_rrtm .EQ. 0) THEN !--old radiation
4037             tau_aero(:,:,:,:) = 1.e-15
4038             piz_aero(:,:,:,:) = 1.
4039             cg_aero(:,:,:,:)  = 0.
4040          ELSE
4041             tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
4042             tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4043             piz_aero_sw_rrtm(:,:,:,:) = 1.0
4044             cg_aero_sw_rrtm(:,:,:,:)  = 0.0
4045          ENDIF
4046       ENDIF
4047       !
[2994]4048       !--WMO criterion to determine tropopause
[3123]4049       CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
[2994]4050       !
[2618]4051       !--STRAT AEROSOL
4052       !--updates tausum_aero,tau_aero,piz_aero,cg_aero
4053       IF (flag_aerosol_strat.GT.0) THEN
4054          IF (prt_level .GE.10) THEN
4055             PRINT *,'appel a readaerosolstrat', mth_cur
4056          ENDIF
4057          IF (iflag_rrtm.EQ.0) THEN
4058           IF (flag_aerosol_strat.EQ.1) THEN
4059             CALL readaerosolstrato(debut)
4060           ELSE
4061             abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
4062             CALL abort_physic(modname,abort_message,1)
4063           ENDIF
4064          ELSE
[2009]4065#ifdef CPP_RRTM
[2690]4066#ifndef CPP_StratAer
4067          !--prescribed strat aerosols
4068          !--only in the case of non-interactive strat aerosols
[2618]4069            IF (flag_aerosol_strat.EQ.1) THEN
4070             CALL readaerosolstrato1_rrtm(debut)
4071            ELSEIF (flag_aerosol_strat.EQ.2) THEN
[3480]4072             CALL readaerosolstrato2_rrtm(debut, ok_volcan)
[2618]4073            ELSE
4074             abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
4075             CALL abort_physic(modname,abort_message,1)
4076            ENDIF
[2690]4077#endif
[2618]4078#else
4079             abort_message='You should compile with -rrtm if running ' &
4080                  // 'with iflag_rrtm=1'
4081             CALL abort_physic(modname,abort_message,1)
4082#endif
4083          ENDIF
[3567]4084       ELSE
4085          tausum_aero(:,:,id_STRAT_phy) = 0.
[2618]4086       ENDIF
[2690]4087!
4088#ifdef CPP_RRTM
4089#ifdef CPP_StratAer
[2692]4090       !--compute stratospheric mask
[3123]4091       CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
[2690]4092       !--interactive strat aerosols
4093       CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
4094#endif
4095#endif
[2618]4096       !--fin STRAT AEROSOL
4097       !     
4098
4099       ! Calculer les parametres optiques des nuages et quelques
4100       ! parametres pour diagnostiques:
4101       !
4102       IF (aerosol_couple.AND.config_inca=='aero') THEN
4103          mass_solu_aero(:,:)    = ccm(:,:,1)
4104          mass_solu_aero_pi(:,:) = ccm(:,:,2)
[2692]4105       ENDIF
[2618]4106
4107       IF (ok_newmicro) then
[3908]4108! AI          IF (iflag_rrtm.NE.0) THEN
4109          IF (iflag_rrtm.EQ.1) THEN
[2618]4110#ifdef CPP_RRTM
4111             IF (ok_cdnc.AND.NRADLP.NE.3) THEN
[2469]4112             abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
4113                  // 'pour ok_cdnc'
[2618]4114             CALL abort_physic(modname,abort_message,1)
4115             ENDIF
[2009]4116#else
4117
[2618]4118             abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
4119             CALL abort_physic(modname,abort_message,1)
[2009]4120#endif
[2618]4121          ENDIF
[3274]4122          CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
[3999]4123               paprs, pplay, t_seri, cldliq, picefra, cldfra, &
[2618]4124               cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
4125               flwp, fiwp, flwc, fiwc, &
4126               mass_solu_aero, mass_solu_aero_pi, &
[4114]4127               cldtaupi, latitude_deg, re, fl, ref_liq, ref_ice, &
[2618]4128               ref_liq_pi, ref_ice_pi)
4129       ELSE
4130          CALL nuage (paprs, pplay, &
[3999]4131               t_seri, cldliq, picefra, cldfra, cldtau, cldemi, &
[2618]4132               cldh, cldl, cldm, cldt, cldq, &
4133               ok_aie, &
4134               mass_solu_aero, mass_solu_aero_pi, &
4135               bl95_b0, bl95_b1, &
4136               cldtaupi, re, fl)
[2469]4137       ENDIF
4138       !
[2618]4139       !IM betaCRF
[2469]4140       !
[2618]4141       cldtaurad   = cldtau
4142       cldtaupirad = cldtaupi
4143       cldemirad   = cldemi
4144       cldfrarad   = cldfra
4145
[2469]4146       !
[2618]4147       IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
4148           lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
4149          !
4150          ! global
4151          !
[3048]4152!IM 251017 begin
[3317]4153!               print*,'physiq betaCRF global zdtime=',zdtime
[3048]4154!IM 251017 end
[2618]4155          DO k=1, klev
4156             DO i=1, klon
4157                IF (pplay(i,k).GE.pfree) THEN
[2469]4158                   beta(i,k) = beta_pbl
[2618]4159                ELSE
[2469]4160                   beta(i,k) = beta_free
[2618]4161                ENDIF
4162                IF (mskocean_beta) THEN
[2469]4163                   beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
[2618]4164                ENDIF
[2469]4165                cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
4166                cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4167                cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
4168                cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
[2618]4169             ENDDO
4170          ENDDO
4171          !
4172       ELSE
4173          !
4174          ! regional
4175          !
4176          DO k=1, klev
4177             DO i=1,klon
4178                !
4179                IF (longitude_deg(i).ge.lon1_beta.AND. &
4180                    longitude_deg(i).le.lon2_beta.AND. &
4181                    latitude_deg(i).le.lat1_beta.AND.  &
4182                    latitude_deg(i).ge.lat2_beta) THEN
4183                   IF (pplay(i,k).GE.pfree) THEN
4184                      beta(i,k) = beta_pbl
4185                   ELSE
4186                      beta(i,k) = beta_free
4187                   ENDIF
4188                   IF (mskocean_beta) THEN
4189                      beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4190                   ENDIF
4191                   cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
4192                   cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4193                   cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
4194                   cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
4195                ENDIF
[2469]4196             !
[2618]4197             ENDDO
[2469]4198          ENDDO
4199       !
[2618]4200       ENDIF
[766]4201
[2618]4202       !lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
4203       IF (ok_chlorophyll) THEN
[2469]4204          print*,"-- reading chlorophyll"
[2618]4205          CALL readchlorophyll(debut)
4206       ENDIF
[1863]4207
[2524]4208!--if ok_suntime_rrtm we use ancillay data for RSUN
4209!--previous values are therefore overwritten
4210!--this is needed for CMIP6 runs
4211!--and only possible for new radiation scheme
4212       IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
[2525]4213#ifdef CPP_RRTM
[2524]4214         CALL read_rsun_rrtm(debut)
[2525]4215#endif
[2524]4216       ENDIF
4217
[2692]4218       IF (mydebug) THEN
4219          CALL writefield_phy('u_seri',u_seri,nbp_lev)
4220          CALL writefield_phy('v_seri',v_seri,nbp_lev)
4221          CALL writefield_phy('t_seri',t_seri,nbp_lev)
4222          CALL writefield_phy('q_seri',q_seri,nbp_lev)
4223       ENDIF
[2524]4224
[2469]4225       !
4226       !sonia : If Iflag_radia >=2, pertubation of some variables
4227       !input to radiation (DICE)
4228       !
4229       IF (iflag_radia .ge. 2) THEN
4230          zsav_tsol (:) = zxtsol(:)
[2692]4231          CALL perturb_radlwsw(zxtsol,iflag_radia)
[2469]4232       ENDIF
[2328]4233
[2469]4234       IF (aerosol_couple.AND.config_inca=='aero') THEN
[959]4235#ifdef INCA
[2469]4236          CALL radlwsw_inca  &
[3338]4237               (chemistry_couple, kdlon,kflev,dist, rmu0, fract, solaire, &
[2469]4238               paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
[2684]4239               size(wo,3), wo, &
[2469]4240               cldfrarad, cldemirad, cldtaurad, &
4241               heat,heat0,cool,cool0,albpla, &
4242               topsw,toplw,solsw,sollw, &
4243               sollwdown, &
4244               topsw0,toplw0,solsw0,sollw0, &
4245               lwdn0, lwdn, lwup0, lwup,  &
4246               swdn0, swdn, swup0, swup, &
4247               ok_ade, ok_aie, &
4248               tau_aero, piz_aero, cg_aero, &
4249               topswad_aero, solswad_aero, &
4250               topswad0_aero, solswad0_aero, &
4251               topsw_aero, topsw0_aero, &
4252               solsw_aero, solsw0_aero, &
4253               cldtaupirad, &
4254               topswai_aero, solswai_aero)
[955]4255#endif
[2469]4256       ELSE
4257          !
4258          !IM calcul radiatif pour le cas actuel
4259          !
4260          RCO2 = RCO2_act
4261          RCH4 = RCH4_act
4262          RN2O = RN2O_act
4263          RCFC11 = RCFC11_act
4264          RCFC12 = RCFC12_act
[3450]4265          !
4266          !--interactive CO2 in ppm from carbon cycle
4267          IF (carbon_cycle_rad.AND..NOT.debut) THEN
4268            RCO2=RCO2_glo
4269          ENDIF
[2469]4270          !
4271          IF (prt_level .GE.10) THEN
4272             print *,' ->radlwsw, number 1 '
4273          ENDIF
4274          !
4275          CALL radlwsw &
4276               (dist, rmu0, fract,  &
4277                                !albedo SB >>>
4278                                !      paprs, pplay,zxtsol,albsol1, albsol2,  &
4279               paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif,  &
4280                                !albedo SB <<<
4281               t_seri,q_seri,wo, &
4282               cldfrarad, cldemirad, cldtaurad, &
[3989]4283               ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, flag_volc_surfstrat, &
4284               flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
[2469]4285               tau_aero, piz_aero, cg_aero, &
4286               tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4287               ! Rajoute par OB pour RRTM
4288               tau_aero_lw_rrtm, &
[3630]4289               cldtaupirad, &
[3048]4290!              zqsat, flwcrad, fiwcrad, &
[2469]4291               zqsat, flwc, fiwc, &
4292               ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4293               heat,heat0,cool,cool0,albpla, &
[3479]4294               heat_volc,cool_volc, &
[3756]4295               topsw,toplw,solsw,solswfdiff,sollw, &
[2469]4296               sollwdown, &
4297               topsw0,toplw0,solsw0,sollw0, &
[3106]4298               lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup,  &
[3082]4299               swdnc0, swdn0, swdn, swupc0, swup0, swup, &
[2469]4300               topswad_aero, solswad_aero, &
4301               topswai_aero, solswai_aero, &
4302               topswad0_aero, solswad0_aero, &
4303               topsw_aero, topsw0_aero, &
4304               solsw_aero, solsw0_aero, &
4305               topswcf_aero, solswcf_aero, &
4306                                !-C. Kleinschmitt for LW diagnostics
4307               toplwad_aero, sollwad_aero,&
4308               toplwai_aero, sollwai_aero, &
4309               toplwad0_aero, sollwad0_aero,&
4310                                !-end
4311               ZLWFT0_i, ZFLDN0, ZFLUP0, &
[3117]4312               ZSWFT0_i, ZFSDN0, ZFSUP0)
[879]4313
[3048]4314          !lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
4315          !schemes
4316          toplw = toplw + betalwoff * (toplw0 - toplw)
4317          sollw = sollw + betalwoff * (sollw0 - sollw)
4318          lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
4319          lwup = lwup + betalwoff * (lwup0 - lwup)
4320          sollwdown(:)= sollwdown(:) + betalwoff *(-1.*ZFLDN0(:,1) - &
4321                        sollwdown(:))
4322          cool = cool + betalwoff * (cool0 - cool)
4323 
[2679]4324#ifndef CPP_XIOS
[4056]4325          !--OB 30/05/2016 modified 21/10/2016
4326          !--here we return swaero_diag and dryaod_diag to FALSE
4327          !--and histdef will switch it back to TRUE if necessary
4328          !--this is necessary to get the right swaero at first step
4329          !--but only in the case of no XIOS as XIOS is covered elsewhere
4330          IF (debut) swaerofree_diag = .FALSE.
4331          IF (debut) swaero_diag = .FALSE.
4332          IF (debut) dryaod_diag = .FALSE.
4333          !--IM 15/09/2017 here we return ok_4xCO2atm to FALSE
4334          !--as for swaero_diag, see above
4335          IF (debut) ok_4xCO2atm = .FALSE.
[3149]4336
[4056]4337          !
[2469]4338          !IM 2eme calcul radiatif pour le cas perturbe ou au moins un
4339          !IM des taux doit etre different du taux actuel
4340          !IM Par defaut on a les taux perturbes egaux aux taux actuels
4341          !
[2989]4342          IF (RCO2_per.NE.RCO2_act.OR. &
4343              RCH4_per.NE.RCH4_act.OR. &
4344              RN2O_per.NE.RN2O_act.OR. &
4345              RCFC11_per.NE.RCFC11_act.OR. &
4346              RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
[3149]4347#endif
[2989]4348   !
[2692]4349          IF (ok_4xCO2atm) THEN
[2469]4350                !
4351                RCO2 = RCO2_per
4352                RCH4 = RCH4_per
4353                RN2O = RN2O_per
4354                RCFC11 = RCFC11_per
4355                RCFC12 = RCFC12_per
4356                !
4357                IF (prt_level .GE.10) THEN
4358                   print *,' ->radlwsw, number 2 '
4359                ENDIF
4360                !
4361                CALL radlwsw &
4362                     (dist, rmu0, fract,  &
4363                                !albedo SB >>>
4364                                !      paprs, pplay,zxtsol,albsol1, albsol2,  &
4365                     paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4366                                !albedo SB <<<
4367                     t_seri,q_seri,wo, &
[2640]4368                     cldfrarad, cldemirad, cldtaurad, &
[3989]4369                     ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, flag_volc_surfstrat, &
4370                     flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
[2469]4371                     tau_aero, piz_aero, cg_aero, &
4372                     tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4373                                ! Rajoute par OB pour RRTM
4374                     tau_aero_lw_rrtm, &
[3630]4375                     cldtaupi, &
[3048]4376!                    zqsat, flwcrad, fiwcrad, &
[2469]4377                     zqsat, flwc, fiwc, &
4378                     ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4379                     heatp,heat0p,coolp,cool0p,albplap, &
[3479]4380                     heat_volc,cool_volc, &
[3756]4381                     topswp,toplwp,solswp,solswfdiffp,sollwp, &
[2469]4382                     sollwdownp, &
4383                     topsw0p,toplw0p,solsw0p,sollw0p, &
[3106]4384                     lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp,  &
[3082]4385                     swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
[2469]4386                     topswad_aerop, solswad_aerop, &
4387                     topswai_aerop, solswai_aerop, &
4388                     topswad0_aerop, solswad0_aerop, &
4389                     topsw_aerop, topsw0_aerop, &
4390                     solsw_aerop, solsw0_aerop, &
4391                     topswcf_aerop, solswcf_aerop, &
4392                                !-C. Kleinschmitt for LW diagnostics
4393                     toplwad_aerop, sollwad_aerop,&
4394                     toplwai_aerop, sollwai_aerop, &
4395                     toplwad0_aerop, sollwad0_aerop,&
4396                                !-end
4397                     ZLWFT0_i, ZFLDN0, ZFLUP0, &
[3117]4398                     ZSWFT0_i, ZFSDN0, ZFSUP0)
[3987]4399          ENDIF !ok_4xCO2atm
[2469]4400       ENDIF ! aerosol_couple
4401       itaprad = 0
4402       !
4403       !  If Iflag_radia >=2, reset pertubed variables
4404       !
4405       IF (iflag_radia .ge. 2) THEN
4406          zxtsol(:) = zsav_tsol (:)
4407       ENDIF
4408    ENDIF ! MOD(itaprad,radpas)
4409    itaprad = itaprad + 1
[879]4410
[2469]4411    IF (iflag_radia.eq.0) THEN
4412       IF (prt_level.ge.9) THEN
4413          PRINT *,'--------------------------------------------------'
4414          PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
4415          PRINT *,'>>>>           heat et cool mis a zero '
4416          PRINT *,'--------------------------------------------------'
[2692]4417       ENDIF
[2469]4418       heat=0.
4419       cool=0.
4420       sollw=0.   ! MPL 01032011
4421       solsw=0.
4422       radsol=0.
4423       swup=0.    ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
4424       swup0=0.
4425       lwup=0.
4426       lwup0=0.
4427       lwdn=0.
4428       lwdn0=0.
[2692]4429    ENDIF
[782]4430
[2469]4431    !
4432    ! Calculer radsol a l'exterieur de radlwsw
4433    ! pour prendre en compte le cycle diurne
4434    ! recode par Olivier Boucher en sept 2015
4435    !
4436    radsol=solsw*swradcorr+sollw
[2618]4437
[2692]4438    IF (ok_4xCO2atm) THEN
[2469]4439       radsolp=solswp*swradcorr+sollwp
[2692]4440    ENDIF
[2359]4441
[2469]4442    !
4443    ! Ajouter la tendance des rayonnements (tous les pas)
4444    ! avec une correction pour le cycle diurne dans le SW
4445    !
[2359]4446
[2469]4447    DO k=1, klev
[3435]4448       d_t_swr(:,k)=swradcorr(:)*heat(:,k)*phys_tstep/RDAY
4449       d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)*phys_tstep/RDAY
4450       d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
4451       d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
[2469]4452    ENDDO
[2194]4453
[2812]4454    CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
[3461]4455    CALL prt_enerbil('SW',itap)
[2812]4456    CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
[3461]4457    CALL prt_enerbil('LW',itap)
[1863]4458
[2469]4459    !
[2692]4460    IF (mydebug) THEN
4461       CALL writefield_phy('u_seri',u_seri,nbp_lev)
4462       CALL writefield_phy('v_seri',v_seri,nbp_lev)
4463       CALL writefield_phy('t_seri',t_seri,nbp_lev)
4464       CALL writefield_phy('q_seri',q_seri,nbp_lev)
4465    ENDIF
[1863]4466
[2469]4467    ! Calculer l'hydrologie de la surface
4468    !
4469    !      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
4470    !     .            agesno, ftsol,fqsurf,fsnow, ruis)
4471    !
[1001]4472
[2469]4473    !
4474    ! Calculer le bilan du sol et la derive de temperature (couplage)
4475    !
4476    DO i = 1, klon
4477       !         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
4478       ! a la demande de JLD
4479       bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
4480    ENDDO
4481    !
4482    !moddeblott(jan95)
4483    ! Appeler le programme de parametrisation de l'orographie
4484    ! a l'echelle sous-maille:
4485    !
4486    IF (prt_level .GE.10) THEN
4487       print *,' call orography ? ', ok_orodr
4488    ENDIF
4489    !
4490    IF (ok_orodr) THEN
4491       !
4492       !  selection des points pour lesquels le shema est actif:
4493       igwd=0
4494       DO i=1,klon
4495          itest(i)=0
4496          !        IF ((zstd(i).gt.10.0)) THEN
4497          IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4498             itest(i)=1
4499             igwd=igwd+1
4500             idx(igwd)=i
4501          ENDIF
4502       ENDDO
4503       !        igwdim=MAX(1,igwd)
4504       !
4505       IF (ok_strato) THEN
[1863]4506
[3435]4507          CALL drag_noro_strato(0,klon,klev,phys_tstep,paprs,pplay, &
[2469]4508               zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4509               igwd,idx,itest, &
4510               t_seri, u_seri, v_seri, &
4511               zulow, zvlow, zustrdr, zvstrdr, &
4512               d_t_oro, d_u_oro, d_v_oro)
[1863]4513
[2469]4514       ELSE
[3435]4515          CALL drag_noro(klon,klev,phys_tstep,paprs,pplay, &
[2469]4516               zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4517               igwd,idx,itest, &
4518               t_seri, u_seri, v_seri, &
4519               zulow, zvlow, zustrdr, zvstrdr, &
4520               d_t_oro, d_u_oro, d_v_oro)
4521       ENDIF
4522       !
4523       !  ajout des tendances
4524       !-----------------------------------------------------------------------
4525       ! ajout des tendances de la trainee de l'orographie
4526       CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
[2812]4527            abortphy,flag_inhib_tend,itap,0)
[3461]4528       CALL prt_enerbil('oro',itap)
[2469]4529       !----------------------------------------------------------------------
4530       !
4531    ENDIF ! fin de test sur ok_orodr
4532    !
[2692]4533    IF (mydebug) THEN
4534       CALL writefield_phy('u_seri',u_seri,nbp_lev)
4535       CALL writefield_phy('v_seri',v_seri,nbp_lev)
4536       CALL writefield_phy('t_seri',t_seri,nbp_lev)
4537       CALL writefield_phy('q_seri',q_seri,nbp_lev)
4538    ENDIF
[1001]4539
[2469]4540    IF (ok_orolf) THEN
4541       !
4542       !  selection des points pour lesquels le shema est actif:
4543       igwd=0
4544       DO i=1,klon
4545          itest(i)=0
4546          IF ((zpic(i)-zmea(i)).GT.100.) THEN
4547             itest(i)=1
4548             igwd=igwd+1
4549             idx(igwd)=i
4550          ENDIF
4551       ENDDO
4552       !        igwdim=MAX(1,igwd)
4553       !
4554       IF (ok_strato) THEN
[1001]4555
[3435]4556          CALL lift_noro_strato(klon,klev,phys_tstep,paprs,pplay, &
[2469]4557               latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
4558               igwd,idx,itest, &
4559               t_seri, u_seri, v_seri, &
4560               zulow, zvlow, zustrli, zvstrli, &
4561               d_t_lif, d_u_lif, d_v_lif               )
[2333]4562
[2469]4563       ELSE
[3435]4564          CALL lift_noro(klon,klev,phys_tstep,paprs,pplay, &
[2469]4565               latitude_deg,zmea,zstd,zpic, &
4566               itest, &
4567               t_seri, u_seri, v_seri, &
4568               zulow, zvlow, zustrli, zvstrli, &
4569               d_t_lif, d_u_lif, d_v_lif)
4570       ENDIF
[1638]4571
[2469]4572       ! ajout des tendances de la portance de l'orographie
4573       CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
[2812]4574            'lif', abortphy,flag_inhib_tend,itap,0)
[3461]4575       CALL prt_enerbil('lif',itap)
[2469]4576    ENDIF ! fin de test sur ok_orolf
[1638]4577
[2469]4578    IF (ok_hines) then
4579       !  HINES GWD PARAMETRIZATION
4580       east_gwstress=0.
4581       west_gwstress=0.
4582       du_gwd_hines=0.
4583       dv_gwd_hines=0.
[3435]4584       CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, &
[2469]4585            u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4586            du_gwd_hines, dv_gwd_hines)
4587       zustr_gwd_hines=0.
4588       zvstr_gwd_hines=0.
4589       DO k = 1, klev
[3435]4590          zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/phys_tstep &
[2469]4591               * (paprs(:, k)-paprs(:, k+1))/rg
[3435]4592          zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/phys_tstep &
[2469]4593               * (paprs(:, k)-paprs(:, k+1))/rg
4594       ENDDO
[1001]4595
[2469]4596       d_t_hin(:, :)=0.
4597       CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
[2812]4598            dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
[3461]4599       CALL prt_enerbil('hin',itap)
[2469]4600    ENDIF
[2333]4601
[2469]4602    IF (.not. ok_hines .and. ok_gwd_rando) then
[3435]4603       ! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod
4604       CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
[2469]4605            v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4606            dv_gwd_front, east_gwstress, west_gwstress)
4607       zustr_gwd_front=0.
4608       zvstr_gwd_front=0.
4609       DO k = 1, klev
[3435]4610          zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/phys_tstep &
[2469]4611               * (paprs(:, k)-paprs(:, k+1))/rg
[3435]4612          zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/phys_tstep &
[2469]4613               * (paprs(:, k)-paprs(:, k+1))/rg
4614       ENDDO
[644]4615
[2469]4616       CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
[2812]4617            paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
[3461]4618       CALL prt_enerbil('front_gwd_rando',itap)
[2469]4619    ENDIF
[1938]4620
[2692]4621    IF (ok_gwd_rando) THEN
[3435]4622       CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
[2469]4623            rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4624            du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4625       CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
[2812]4626            paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
[3461]4627       CALL prt_enerbil('flott_gwd_rando',itap)
[2469]4628       zustr_gwd_rando=0.
4629       zvstr_gwd_rando=0.
4630       DO k = 1, klev
[3435]4631          zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/phys_tstep &
[2469]4632               * (paprs(:, k)-paprs(:, k+1))/rg
[3435]4633          zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/phys_tstep &
[2469]4634               * (paprs(:, k)-paprs(:, k+1))/rg
4635       ENDDO
[2692]4636    ENDIF
[766]4637
[2469]4638    ! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
[1279]4639
[2692]4640    IF (mydebug) THEN
4641       CALL writefield_phy('u_seri',u_seri,nbp_lev)
4642       CALL writefield_phy('v_seri',v_seri,nbp_lev)
4643       CALL writefield_phy('t_seri',t_seri,nbp_lev)
4644       CALL writefield_phy('q_seri',q_seri,nbp_lev)
4645    ENDIF
[2136]4646
[2469]4647    DO i = 1, klon
4648       zustrph(i)=0.
4649       zvstrph(i)=0.
4650    ENDDO
4651    DO k = 1, klev
4652       DO i = 1, klon
[3435]4653          zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/phys_tstep* &
[2469]4654               (paprs(i,k)-paprs(i,k+1))/rg
[3435]4655          zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/phys_tstep* &
[2469]4656               (paprs(i,k)-paprs(i,k+1))/rg
4657       ENDDO
4658    ENDDO
4659    !
4660    !IM calcul composantes axiales du moment angulaire et couple des montagnes
4661    !
4662    IF (is_sequential .and. ok_orodr) THEN
4663       CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
4664            ra,rg,romega, &
4665            latitude_deg,longitude_deg,pphis, &
4666            zustrdr,zustrli,zustrph, &
4667            zvstrdr,zvstrli,zvstrph, &
4668            paprs,u,v, &
4669            aam, torsfc)
4670    ENDIF
4671    !IM cf. FLott END
4672    !DC Calcul de la tendance due au methane
[3461]4673    IF (ok_qch4) THEN
[2469]4674       CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
4675       ! ajout de la tendance d'humidite due au methane
[3435]4676       d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
[2801]4677       CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
[2812]4678            'q_ch4', abortphy,flag_inhib_tend,itap,0)
[3435]4679       d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep
[2692]4680    ENDIF
[2469]4681    !
4682    !
[2897]4683
4684!===============================================================
4685!            Additional tendency of TKE due to orography
4686!===============================================================
4687!
4688! Inititialization
4689!------------------
4690
4691       addtkeoro=0   
4692       CALL getin_p('addtkeoro',addtkeoro)
4693     
4694       IF (prt_level.ge.5) &
4695            print*,'addtkeoro', addtkeoro
4696           
4697       alphatkeoro=1.   
4698       CALL getin_p('alphatkeoro',alphatkeoro)
4699       alphatkeoro=min(max(0.,alphatkeoro),1.)
4700
[3461]4701       smallscales_tkeoro=.FALSE.   
[2897]4702       CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
4703
4704
[3461]4705       dtadd(:,:)=0.
4706       duadd(:,:)=0.
4707       dvadd(:,:)=0.
[2897]4708
4709! Choices for addtkeoro:
4710!      ** 0 no TKE tendency from orography   
4711!      ** 1 we include a fraction alphatkeoro of the whole tendency duoro
4712!      ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
4713!
4714
4715       IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
4716!      -------------------------------------------
4717
4718
4719       !  selection des points pour lesquels le schema est actif:
4720
4721
4722  IF (addtkeoro .EQ. 1 ) THEN
4723
4724            duadd(:,:)=alphatkeoro*d_u_oro(:,:)
4725            dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
4726
4727  ELSE IF (addtkeoro .EQ. 2) THEN
4728
[3461]4729     IF (smallscales_tkeoro) THEN
[2897]4730       igwd=0
4731       DO i=1,klon
4732          itest(i)=0
4733! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
4734! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
4735! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
4736          IF (zstd(i).GT.1.0) THEN
4737             itest(i)=1
4738             igwd=igwd+1
4739             idx(igwd)=i
4740          ENDIF
4741       ENDDO
4742
4743     ELSE
4744
4745       igwd=0
4746       DO i=1,klon
4747          itest(i)=0
4748        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4749             itest(i)=1
4750             igwd=igwd+1
4751             idx(igwd)=i
[3461]4752        ENDIF
[2897]4753       ENDDO
4754
[3461]4755     ENDIF
[2897]4756
[3461]4757     CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
[2897]4758               zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4759               igwd,idx,itest, &
4760               t_seri, u_seri, v_seri, &
4761               zulow, zvlow, zustrdr, zvstrdr, &
4762               d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
4763
[3461]4764     zustrdr(:)=0.
4765     zvstrdr(:)=0.
4766     zulow(:)=0.
4767     zvlow(:)=0.
[2897]4768
[3461]4769     duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
4770     dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
4771  ENDIF
[2897]4772
4773
4774   ! TKE update from subgrid temperature and wind tendencies
4775   !----------------------------------------------------------
4776    forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4777
4778
[3198]4779    CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
[3888]4780   !
4781   ! Prevent pbl_tke_w from becoming negative
4782    wake_delta_pbl_tke(:,:,:) = max(wake_delta_pbl_tke(:,:,:), -pbl_tke(:,:,:))
4783   !
[2897]4784
4785       ENDIF
4786!      -----
4787!===============================================================
4788
4789
[2469]4790    !====================================================================
4791    ! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4792    !====================================================================
4793    ! Abderrahmane 24.08.09
4794
4795    IF (ok_cosp) THEN
4796       ! adeclarer
[1279]4797#ifdef CPP_COSP
[3435]4798       IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
[1279]4799
[2469]4800          IF (prt_level .GE.10) THEN
4801             print*,'freq_cosp',freq_cosp
4802          ENDIF
4803          mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4804          !       print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4805          !     s        ref_liq,ref_ice
[3435]4806          CALL phys_cosp(itap,phys_tstep,freq_cosp, &
[2469]4807               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
[2794]4808               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
[2469]4809               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4810               JrNt,ref_liq,ref_ice, &
4811               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4812               zu10m,zv10m,pphis, &
4813               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4814               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4815               prfl(:,1:klev),psfl(:,1:klev), &
4816               pmflxr(:,1:klev),pmflxs(:,1:klev), &
4817               mr_ozone,cldtau, cldemi)
[1412]4818
[2469]4819          !     L         calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4820          !     L          cfaddbze,clcalipso2,dbze,cltlidarradar,
4821          !     M          clMISR,
4822          !     R          clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4823          !     I          tauisccp,albisccp,meantbisccp,meantbclrisccp)
[1279]4824
[2469]4825       ENDIF
[3370]4826#endif
[1279]4827
[3370]4828#ifdef CPP_COSP2
[3435]4829       IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
[3370]4830
4831          IF (prt_level .GE.10) THEN
4832             print*,'freq_cosp',freq_cosp
4833          ENDIF
4834          mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4835                 print*,'Dans physiq.F avant appel '
4836          !     s        ref_liq,ref_ice
[3435]4837          CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
[3370]4838               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4839               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4840               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4841               JrNt,ref_liq,ref_ice, &
4842               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4843               zu10m,zv10m,pphis, &
4844               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4845               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4846               prfl(:,1:klev),psfl(:,1:klev), &
4847               pmflxr(:,1:klev),pmflxs(:,1:klev), &
4848               mr_ozone,cldtau, cldemi)
4849       ENDIF
[1279]4850#endif
[3370]4851
[3491]4852#ifdef CPP_COSPV2
4853       IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
[3726]4854!        IF (MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
[3491]4855
4856          IF (prt_level .GE.10) THEN
4857             print*,'freq_cosp',freq_cosp
4858          ENDIF
[3726]4859           DO k = 1, klev
4860             DO i = 1, klon
4861               phicosp(i,k) = pphi(i,k) + pphis(i)
4862             ENDDO
4863           ENDDO
[3491]4864          mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4865                 print*,'Dans physiq.F avant appel '
4866          !     s        ref_liq,ref_ice
4867          CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
4868               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4869               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4870               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4871               JrNt,ref_liq,ref_ice, &
4872               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4873               zu10m,zv10m,pphis, &
4874               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4875               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4876               prfl(:,1:klev),psfl(:,1:klev), &
4877               pmflxr(:,1:klev),pmflxs(:,1:klev), &
4878               mr_ozone,cldtau, cldemi)
4879       ENDIF
4880#endif
4881
[2469]4882    ENDIF  !ok_cosp
[2580]4883
4884
4885! Marine
4886
4887  IF (ok_airs) then
4888
[3435]4889  IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
[2692]4890     write(*,*) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
4891     CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
4892        & map_prop_hc,map_prop_hist,&
4893        & map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
4894        & map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
4895        & map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
4896        & map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
4897        & map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
4898        & map_ntot,map_hc,map_hist,&
4899        & map_Cb,map_ThCi,map_Anv,&
4900        & alt_tropo )
[2580]4901  ENDIF
4902
4903  ENDIF  ! ok_airs
4904
4905
[2469]4906    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4907    !AA
4908    !AA Installation de l'interface online-offline pour traceurs
4909    !AA
4910    !====================================================================
4911    !   Calcul  des tendances traceurs
4912    !====================================================================
4913    !
[959]4914
[2469]4915    IF (type_trac=='repr') THEN
[3666]4916!MM pas d'impact, car on recupere q_seri,tr_seri,t_seri via phys_local_var_mod
4917!MM                               dans Reprobus
[2469]4918       sh_in(:,:) = q_seri(:,:)
[3666]4919#ifdef REPROBUS
4920       d_q_rep(:,:) = 0.
4921       d_ql_rep(:,:) = 0.
4922       d_qi_rep(:,:) = 0.
4923#endif
[2469]4924    ELSE
4925       sh_in(:,:) = qx(:,:,ivap)
[4098]4926       IF (nqo >= 3) THEN
[3861]4927          ch_in(:,:) = qx(:,:,iliq) + qx(:,:,isol)
4928       ELSE
4929          ch_in(:,:) = qx(:,:,iliq)
4930       ENDIF
[2692]4931    ENDIF
[1565]4932
[2630]4933#ifdef CPP_Dust
[3776]4934    !  Avec SPLA, iflag_phytrac est forcé =1
4935    CALL       phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con,       &  ! I
[2630]4936                      pdtphys,ftsol,                                   &  ! I
4937                      t,q_seri,paprs,pplay,RHcl,                  &  ! I
4938                      pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,          &  ! I
4939                      coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1,                 &  ! I
4940                      u_seri, v_seri, latitude_deg, longitude_deg,  &
4941                      pphis,pctsrf,pmflxr,pmflxs,prfl,psfl,            &  ! I
4942                      da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij,     &  ! I
4943                      epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con,      &  ! I
4944                      ev,wdtrainA,  wdtrainM,wght_cvfd,              &  ! I
4945                      fm_therm, entr_therm, rneb,                      &  ! I
4946                      beta_prec_fisrt,beta_prec, & !I
4947                      zu10m,zv10m,wstar,ale_bl,ale_wake,               &  ! I
4948                      d_tr_dyn,tr_seri)
4949
4950#else
[3776]4951    IF (iflag_phytrac == 1 ) THEN
4952      CALL phytrac ( &
[2469]4953         itap,     days_elapsed+1,    jH_cur,   debut, &
[3435]4954         lafin,    phys_tstep,     u, v,     t, &
[2469]4955         paprs,    pplay,     pmfu,     pmfd, &
4956         pen_u,    pde_u,     pen_d,    pde_d, &
4957         cdragh,   coefh(1:klon,1:klev,is_ave),   fm_therm, entr_therm, &
4958         u1,       v1,        ftsol,    pctsrf, &
4959         zustar,   zu10m,     zv10m, &
4960         wstar(:,is_ave),    ale_bl,         ale_wake, &
4961         latitude_deg, longitude_deg, &
4962         frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4963         presnivs, pphis,     pphi,     albsol1, &
[2784]4964         sh_in,   ch_in,    rhcl,      cldfra,   rneb, &
[2469]4965         diafra,   cldliq,    itop_con, ibas_con, &
4966         pmflxr,   pmflxs,    prfl,     psfl, &
4967         da,       phi,       mp,       upwd, &
4968         phi2,     d1a,       dam,      sij, wght_cvfd, &        !<<RomP+RL
4969         wdtrainA, wdtrainM,  sigd,     clw,elij, &   !<<RomP
4970         ev,       ep,        epmlmMm,  eplaMm, &     !<<RomP
4971         dnwd,     aerosol_couple,      flxmass_w, &
4972         tau_aero, piz_aero,  cg_aero,  ccm, &
4973         rfname, &
4974         d_tr_dyn, &                                 !<<RomP
[3418]4975         tr_seri, init_source)
[3666]4976#ifdef REPROBUS
4977
4978
4979          print*,'avt add phys rep',abortphy
4980
4981     CALL add_phys_tend &
4982            (du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,paprs,&
4983             'rep',abortphy,flag_inhib_tend,itap,0)
4984        IF (abortphy==1) Print*,'ERROR ABORT REP'
4985
4986          print*,'apr add phys rep',abortphy
4987
[2630]4988#endif
[3776]4989    ENDIF    ! (iflag_phytrac=1)
[3666]4990
4991#endif
[3776]4992    !ENDIF    ! (iflag_phytrac=1)
[524]4993
[2469]4994    IF (offline) THEN
[524]4995
[2469]4996       IF (prt_level.ge.9) &
4997            print*,'Attention on met a 0 les thermiques pour phystoke'
[2692]4998       CALL phystokenc ( &
[2469]4999            nlon,klev,pdtphys,longitude_deg,latitude_deg, &
5000            t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
5001            fm_therm,entr_therm, &
5002            cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
5003            frac_impa, frac_nucl, &
[3435]5004            pphis,cell_area,phys_tstep,itap, &
[2469]5005            qx(:,:,ivap),da,phi,mp,upwd,dnwd)
[524]5006
5007
[2469]5008    ENDIF
[524]5009
[2469]5010    !
5011    ! Calculer le transport de l'eau et de l'energie (diagnostique)
5012    !
5013    CALL transp (paprs,zxtsol, &
[3257]5014         t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
5015         ve, vq, ue, uq, vwat, uwat)
[2469]5016    !
5017    !IM global posePB BEG
5018    IF(1.EQ.0) THEN
5019       !
5020       CALL transp_lay (paprs,zxtsol, &
5021            t_seri, q_seri, u_seri, v_seri, zphi, &
5022            ve_lay, vq_lay, ue_lay, uq_lay)
5023       !
5024    ENDIF !(1.EQ.0) THEN
5025    !IM global posePB END
5026    ! Accumuler les variables a stocker dans les fichiers histoire:
5027    !
[1279]5028
[2469]5029    !================================================================
5030    ! Conversion of kinetic and potential energy into heat, for
5031    ! parameterisation of subgrid-scale motions
5032    !================================================================
[1753]5033
[2469]5034    d_t_ec(:,:)=0.
5035    forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
[2851]5036    CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
5037         u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
[2469]5038         zmasse,exner,d_t_ec)
5039    t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
[1753]5040
[2469]5041    !=======================================================================
5042    !   SORTIES
5043    !=======================================================================
5044    !
5045    !IM initialisation + calculs divers diag AMIP2
5046    !
5047    include "calcul_divers.h"
5048    !
5049    !IM Interpolation sur les niveaux de pression du NMC
5050    !   -------------------------------------------------
5051    !
5052    include "calcul_STDlev.h"
5053    !
5054    ! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
5055    CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
5056    !
[2496]5057    !cc prw  = eau precipitable
5058    !   prlw = colonne eau liquide
5059    !   prlw = colonne eau solide
[2499]5060    prw(:) = 0.
5061    prlw(:) = 0.
5062    prsw(:) = 0.
5063    DO k = 1, klev
5064       prw(:)  = prw(:)  + q_seri(:,k)*zmasse(:,k)
5065       prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
5066       prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
[2469]5067    ENDDO
5068    !
[3872]5069    IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN
[655]5070#ifdef INCA
[2469]5071       CALL VTe(VTphysiq)
5072       CALL VTb(VTinca)
[959]5073
[2469]5074       CALL chemhook_end ( &
[3435]5075            phys_tstep, &
[2469]5076            pplay, &
5077            t_seri, &
[3872]5078            tr_seri(:,:,1+nqCO2:nbtr), &
[2469]5079            nbtr, &
5080            paprs, &
5081            q_seri, &
5082            cell_area, &
5083            pphi, &
5084            pphis, &
[2832]5085            zx_rh, &
[3632]5086            aps, bps, ap, bp)
[959]5087
[2469]5088       CALL VTe(VTinca)
5089       CALL VTb(VTphysiq)
[655]5090#endif
[2692]5091    ENDIF
[655]5092
[1753]5093
[2469]5094    !
5095    ! Convertir les incrementations en tendances
5096    !
5097    IF (prt_level .GE.10) THEN
5098       print *,'Convertir les incrementations en tendances '
5099    ENDIF
5100    !
[2692]5101    IF (mydebug) THEN
5102       CALL writefield_phy('u_seri',u_seri,nbp_lev)
5103       CALL writefield_phy('v_seri',v_seri,nbp_lev)
5104       CALL writefield_phy('t_seri',t_seri,nbp_lev)
5105       CALL writefield_phy('q_seri',q_seri,nbp_lev)
5106    ENDIF
[766]5107
[2469]5108    DO k = 1, klev
5109       DO i = 1, klon
[3435]5110          d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
5111          d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
5112          d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
5113          d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
5114          d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
[2469]5115          !CR: on ajoute le contenu en glace
[4098]5116          IF (nqo >= 3) THEN
[3435]5117             d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
[2692]5118          ENDIF
[4059]5119          !--ice_sursat: nqo=4, on ajoute rneb
[4098]5120          IF (nqo == 4) THEN
[4059]5121             d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
5122          ENDIF
[2469]5123       ENDDO
5124    ENDDO
5125    !
[4069]5126    IF (nqtot > nqo) THEN
[4056]5127       itr = 0
5128       DO iq = 1, nqtot
[4071]5129          IF(.NOT.tracers(iq)%isInPhysics) CYCLE
[4056]5130          itr = itr+1
[2469]5131          DO  k = 1, klev
5132             DO  i = 1, klon
[4056]5133                d_qx(i,k,iq) = ( tr_seri(i,k,itr) - qx(i,k,iq) ) / phys_tstep
[2469]5134             ENDDO
5135          ENDDO
5136       ENDDO
5137    ENDIF
5138    !
5139    !IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
5140    !IM global posePB      include "write_bilKP_ins.h"
5141    !IM global posePB      include "write_bilKP_ave.h"
5142    !
[1412]5143
[2489]5144    !--OB mass fixer
5145    !--profile is corrected to force mass conservation of water
5146    IF (mass_fixer) THEN
5147    qql2(:)=0.0
[2499]5148    DO k = 1, klev
5149      qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
[2489]5150    ENDDO
5151    DO i = 1, klon
5152      !--compute ratio of what q+ql should be with conservation to what it is
5153      corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
5154      DO k = 1, klev
5155        q_seri(i,k) =q_seri(i,k)*corrqql
5156        ql_seri(i,k)=ql_seri(i,k)*corrqql
5157      ENDDO
5158    ENDDO
5159    ENDIF
5160    !--fin mass fixer
5161
[2469]5162    ! Sauvegarder les valeurs de t et q a la fin de la physique:
5163    !
[2499]5164    u_ancien(:,:)  = u_seri(:,:)
5165    v_ancien(:,:)  = v_seri(:,:)
5166    t_ancien(:,:)  = t_seri(:,:)
5167    q_ancien(:,:)  = q_seri(:,:)
5168    ql_ancien(:,:) = ql_seri(:,:)
5169    qs_ancien(:,:) = qs_seri(:,:)
[4059]5170    rneb_ancien(:,:) = rneb_seri(:,:)
[2499]5171    CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
5172    CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
5173    CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
[2469]5174    ! !! RomP >>>
[4056]5175    IF (nqtot > nqo) tr_ancien(:,:,:) = tr_seri(:,:,:)
[2469]5176    ! !! RomP <<<
5177    !==========================================================================
5178    ! Sorties des tendances pour un point particulier
5179    ! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
5180    ! pour le debug
5181    ! La valeur de igout est attribuee plus haut dans le programme
5182    !==========================================================================
[879]5183
[2692]5184    IF (prt_level.ge.1) THEN
[2469]5185       write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
5186       write(lunout,*) &
5187            'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
5188       write(lunout,*) &
5189            nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
5190            pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
5191            pctsrf(igout,is_sic)
5192       write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
[2692]5193       DO k=1,klev
[2469]5194          write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
5195               d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
5196               d_t_eva(igout,k)
[2692]5197       ENDDO
[2469]5198       write(lunout,*) 'cool,heat'
[2692]5199       DO k=1,klev
[2469]5200          write(lunout,*) cool(igout,k),heat(igout,k)
[2692]5201       ENDDO
[879]5202
[2469]5203       !jyg<     (En attendant de statuer sur le sort de d_t_oli)
5204       !jyg!     write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5205       !jyg!     do k=1,klev
5206       !jyg!        write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
5207       !jyg!             d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5208       !jyg!     enddo
5209       write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
[2692]5210       DO k=1,klev
[2469]5211          write(lunout,*) d_t_vdf(igout,k), &
5212               d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
[2692]5213       ENDDO
[2469]5214       !>jyg
[879]5215
[2469]5216       write(lunout,*) 'd_ps ',d_ps(igout)
5217       write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
[2692]5218       DO k=1,klev
[2469]5219          write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
5220               d_qx(igout,k,1),d_qx(igout,k,2)
[2692]5221       ENDDO
5222    ENDIF
[879]5223
[2469]5224    !============================================================
5225    !   Calcul de la temperature potentielle
5226    !============================================================
5227    DO k = 1, klev
5228       DO i = 1, klon
5229          !JYG/IM theta en debut du pas de temps
5230          !JYG/IM       theta(i,k)=t(i,k)*(100000./pplay(i,k))**(RD/RCPD)
5231          !JYG/IM theta en fin de pas de temps de physique
5232          theta(i,k)=t_seri(i,k)*(100000./pplay(i,k))**(RD/RCPD)
5233          ! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
5234          !     MPL 20130625
5235          ! fth_fonctions.F90 et parkind1.F90
5236          ! sinon thetal=theta
5237          !       thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
5238          !    :         ql_seri(i,k))
5239          thetal(i,k)=theta(i,k)
5240       ENDDO
5241    ENDDO
5242    !
[879]5243
[2469]5244    ! 22.03.04 BEG
5245    !=============================================================
5246    !   Ecriture des sorties
5247    !=============================================================
[524]5248#ifdef CPP_IOIPSL
5249
[2469]5250    ! Recupere des varibles calcule dans differents modules
5251    ! pour ecriture dans histxxx.nc
[782]5252
[2469]5253    ! Get some variables from module fonte_neige_mod
5254    CALL fonte_neige_get_vars(pctsrf,  &
[2517]5255         zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
[782]5256
[1507]5257
[2469]5258    !=============================================================
5259    ! Separation entre thermiques et non thermiques dans les sorties
5260    ! de fisrtilp
5261    !=============================================================
[1507]5262
[2692]5263    IF (iflag_thermals>=1) THEN
[2469]5264       d_t_lscth=0.
5265       d_t_lscst=0.
5266       d_q_lscth=0.
5267       d_q_lscst=0.
[2692]5268       DO k=1,klev
5269          DO i=1,klon
5270             IF (ptconvth(i,k)) THEN
[2469]5271                d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5272                d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
[2692]5273             ELSE
[2469]5274                d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5275                d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
[2692]5276             ENDIF
5277          ENDDO
5278       ENDDO
[1507]5279
[2692]5280       DO i=1,klon
[2469]5281          plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
5282          plul_th(i)=prfl(i,1)+psfl(i,1)
[2692]5283       ENDDO
5284    ENDIF
[909]5285
[2469]5286    !On effectue les sorties:
[1791]5287
[2630]5288#ifdef CPP_Dust
5289  CALL phys_output_write_spl(itap, pdtphys, paprs, pphis,  &
5290       pplay, lmax_th, aerosol_couple,                 &
[3776]5291       ok_ade, ok_aie, ivap, ok_sync,                  &
[2630]5292       ptconv, read_climoz, clevSTD,                   &
5293       ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse,      &
5294       flag_aerosol, flag_aerosol_strat, ok_cdnc)
5295#else
[2469]5296    CALL phys_output_write(itap, pdtphys, paprs, pphis,  &
5297         pplay, lmax_th, aerosol_couple,                 &
[3630]5298         ok_ade, ok_aie, ok_volcan, ivap, iliq, isol,    &
[2496]5299         ok_sync, ptconv, read_climoz, clevSTD,          &
[2665]5300         ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
[2469]5301         flag_aerosol, flag_aerosol_strat, ok_cdnc)
[2630]5302#endif
[1791]5303
[2651]5304#ifndef CPP_XIOS
[2590]5305    CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
[2651]5306#endif
[687]5307
[524]5308#endif
5309
[2469]5310    !====================================================================
5311    ! Arret du modele apres hgardfou en cas de detection d'un
5312    ! plantage par hgardfou
5313    !====================================================================
[2235]5314
5315    IF (abortphy==1) THEN
5316       abort_message ='Plantage hgardfou'
[2311]5317       CALL abort_physic (modname,abort_message,1)
[2235]5318    ENDIF
5319
[2469]5320    ! 22.03.04 END
5321    !
5322    !====================================================================
5323    ! Si c'est la fin, il faut conserver l'etat de redemarrage
5324    !====================================================================
5325    !
[782]5326
[3981]5327    ! Disabling calls to the prt_alerte function
5328    alert_first_call = .FALSE.
5329   
[2469]5330    IF (lafin) THEN
5331       itau_phy = itau_phy + itap
5332       CALL phyredem ("restartphy.nc")
5333       !         open(97,form="unformatted",file="finbin")
5334       !         write(97) u_seri,v_seri,t_seri,q_seri
5335       !         close(97)
[3435]5336     
5337       IF (is_omp_master) THEN
5338       
5339         IF (read_climoz >= 1) THEN
5340           IF (is_mpi_root) CALL nf95_close(ncid_climoz)
5341            DEALLOCATE(press_edg_climoz) ! pointer
5342            DEALLOCATE(press_cen_climoz) ! pointer
5343         ENDIF
5344       
[2692]5345       ENDIF
[3435]5346#ifdef CPP_XIOS
5347       IF (is_omp_master) CALL xios_context_finalize
5348#endif
[3461]5349       WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
[2469]5350    ENDIF
[1863]5351
[2469]5352    !      first=.false.
[1863]5353
[2469]5354  END SUBROUTINE physiq
[2418]5355
[2902]5356END MODULE physiq_mod
Note: See TracBrowser for help on using the repository browser.