source: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90 @ 3256

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

XIOS/OMP problem resolved

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