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

Last change on this file since 2660 was 2657, checked in by fhourdin, 8 years ago

Seuil max sur eau glace nuageuse oicemax (en plus de oliqmax)
Truc pour éviter des plantages dans des atmosphère super saturées.

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