source: LMDZ6/branches/SETHET_DECOUPLE/libf/phylmd/physiq_mod.F90 @ 5182

Last change on this file since 5182 was 3823, checked in by musat, 4 years ago

Nouveaux calculs a 2m et 10m

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