source: LMDZ5/branches/testing/libf/phylmd/physiq_mod.F90 @ 5212

Last change on this file since 5212 was 2886, checked in by Laurent Fairhead, 7 years ago

Merged trunk changes r2865:2885 into testing branch

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