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

Last change on this file since 2661 was 2661, checked in by Laurent Fairhead, 8 years ago

Bug correction on the ozone reading:
When read_climoz = 1 or 2, values for ozone are read from a file with 360 daily values
so that a special calendar is used to read in the proper values when using a 365/366
day calendar in the model. Problem was the same special calendar was used to define
when the ozone was supposed to be updated resulting in an update that became
desynchronised with the actual change of days in the simulation. (is this making sense?)
Took advantage of this correction to rewrite the ozone update mechanism
LF

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