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

Last change on this file since 2606 was 2606, checked in by jyg, 8 years ago

Change of arguments of 'calltherm' in 'physiq' so
that t_seri, q_seri, ... variables are changed by
a call to 'add_phys_tend'.

The 'USE phys_local_var_mod' instruction in
'physiq' is completed with 'ONLY' and the list of
the relevant variables.

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