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

Last change on this file since 3409 was 3387, checked in by oboucher, 6 years ago

adding infocfields_init routine and
its call from physiq. This routines reads
which fields need to be transferred between
model components in ESM configuration.
No impact whatsoever in LMDZ mode.

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