source: LMDZ5/trunk/libf/phylmd/physiq_mod.F90 @ 2648

Last change on this file since 2648 was 2644, checked in by oboucher, 8 years ago

I introduced flag_bc_internal_mixture for BC/sulphate internal mixture.
Only works for iflag_rrtm=y, NSW=6 and flag_aerosol=6 or aerosol_couple.
It has no impact at the moment as LUT for aerosol properties will have to
be changed. But the plumbing work is already done....

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