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

Last change on this file since 2618 was 2618, checked in by oboucher, 8 years ago

I have moved the computation of cloud optical properties
into the radiation if .. .endif sequence, after the aerosol
optical properties. This follows up on rev 2614 and should
fix the 1+1=2 problem that was introduced back then.

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