source: LMDZ5/branches/IPSLCM6.0.8/libf/phylmd/physiq_mod.F90 @ 2791

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

Merged trunk changes r2727:2785 into testing branch

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