source: LMDZ5/branches/testing/libf/phylmd/physiq.F90 @ 2187

Last change on this file since 2187 was 2187, checked in by Laurent Fairhead, 9 years ago

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