source: LMDZ5/branches/IPSLCM6.0.11/libf/phylmd/physiq_mod.F90 @ 3717

Last change on this file since 3717 was 2906, checked in by acozic, 7 years ago

Add an argument to initialisation of INCA

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