source: LMDZ5/branches/IPSLCM6.0.10/libf/phylmd/physiq_mod.F90 @ 5445

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

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