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

Last change on this file since 3600 was 3600, checked in by oboucher, 5 years ago

adding a ompthreadprivate to ok_volcan

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