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

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

Printing the current time and date when entering the physics (information is needed when
extracting crashing points). Date is printed out by mpi/omp master process every 5 physics
steps or at each physics timestep if prt_level > 5

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