source: LMDZ5/trunk/libf/phylmd/physiq.F90 @ 2224

Last change on this file since 2224 was 2224, checked in by crio, 9 years ago

cv3p1_closure: Bug fix in the toothpaste closure computation (dimension of klfc).
physiq: check number of H2O tracers if thermodynamical effect of ice is included. Stop simulation if ice water is missing.

  • 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.8 KB
Line 
1! $Id: physiq.F90 2224 2015-03-11 11:23:40Z crio $
2!#define IO_DEBUG
3
4SUBROUTINE physiq (nlon,nlev, &
5     debut,lafin,jD_cur, jH_cur,pdtphys, &
6     paprs,pplay,pphi,pphis,presnivs, &
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,PARAMETER :: longcles=20
286  REAL,SAVE :: clesphy0(longcles)
287  !$OMP THREADPRIVATE(clesphy0)
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_cldth
639  save iflag_cldth
640  !$OMP THREADPRIVATE(iflag_cldth)
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_cldth,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!CR: check sur le nb de traceurs de l eau
927     if ((iflag_ice_thermo.gt.0).and.(nqo==2)) then
928          WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers (H2Ov, H2Ol, H2Oi)', ' but nqo=', nqo, &
929          '. Might as well stop here.'
930          STOP
931     endif
932
933     dnwd0=0.0
934     ftd=0.0
935     fqd=0.0
936     cin=0.
937     !ym Attention pbase pas initialise dans concvl !!!!
938     pbase=0
939     !IM 180608
940
941     itau_con=0
942     first=.false.
943
944  endif  ! first
945
946  !ym => necessaire pour iflag_con != 2   
947  pmfd(:,:) = 0.
948  pen_u(:,:) = 0.
949  pen_d(:,:) = 0.
950  pde_d(:,:) = 0.
951  pde_u(:,:) = 0.
952  aam=0.
953
954  torsfc=0.
955  forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
956
957
958
959  modname = 'physiq'
960  !IM
961  IF (ip_ebil_phy.ge.1) THEN
962     DO i=1,klon
963        zero_v(i)=0.
964     END DO
965  END IF
966
967  IF (debut) THEN
968     CALL suphel ! initialiser constantes et parametres phys.
969  ENDIF
970
971  if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
972
973
974  !======================================================================
975  ! Gestion calendrier : mise a jour du module phys_cal_mod
976  !
977  !     CALL phys_cal_update(jD_cur,jH_cur)
978
979  !
980  ! Si c'est le debut, il faut initialiser plusieurs choses
981  !          ********
982  !
983  IF (debut) THEN
984     !rv
985     !CRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
986     !de la convection a partir des caracteristiques du thermique
987     wght_th(:,:)=1.
988     lalim_conv(:)=1
989     !RC
990     ustar(:,:)=0.
991     u10m(:,:)=0.
992     v10m(:,:)=0.
993     rain_con(:)=0.
994     snow_con(:)=0.
995     topswai(:)=0.
996     topswad(:)=0.
997     solswai(:)=0.
998     solswad(:)=0.
999
1000     wmax_th(:)=0.
1001     tau_overturning_th(:)=0.
1002
1003     IF (type_trac == 'inca') THEN
1004        ! jg : initialisation jusqu'au ces variables sont dans restart
1005        ccm(:,:,:) = 0.
1006        tau_aero(:,:,:,:) = 0.
1007        piz_aero(:,:,:,:) = 0.
1008        cg_aero(:,:,:,:) = 0.
1009     END IF
1010
1011     rnebcon0(:,:) = 0.0
1012     clwcon0(:,:) = 0.0
1013     rnebcon(:,:) = 0.0
1014     clwcon(:,:) = 0.0
1015
1016     !IM     
1017     IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1018     !
1019     print*,'iflag_coupl,iflag_clos,iflag_wake', &
1020          iflag_coupl,iflag_clos,iflag_wake
1021     print*,'CYCLE_DIURNE', cycle_diurne
1022     !
1023     IF (iflag_con.EQ.2.AND.iflag_cldth.GT.-1) THEN
1024        abort_message = 'Tiedtke needs iflag_cldth=-2 or -1'
1025        CALL abort_gcm (modname,abort_message,1)
1026     ENDIF
1027     !
1028     !
1029     ! Initialiser les compteurs:
1030     !
1031     itap    = 0
1032     itaprad = 0
1033
1034!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1035     !! Un petit travail \`a faire ici.
1036!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1037
1038     if (iflag_pbl>1) then
1039        PRINT*, "Using method MELLOR&YAMADA"
1040     endif
1041
1042!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1043     ! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1044     ! dyn3d
1045     ! Attention : la version precedente n'etait pas tres propre.
1046     ! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1047     ! pour obtenir le meme resultat.
1048     dtime=pdtphys
1049     radpas = NINT( 86400./dtime/nbapp_rad)
1050!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1051
1052     CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1053     IF (klon_glo==1) THEN
1054        coefh=0. ; coefm=0. ; pbl_tke=0.
1055        coefh(:,2,:)=1.e-2 ; coefm(:,2,:)=1.e-2 ; pbl_tke(:,2,:)=1.e-2
1056        PRINT*,'FH WARNING : lignes a supprimer'
1057     ENDIF
1058     !IM begin
1059     print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1060          ,ratqs(1,1)
1061     !IM end
1062
1063
1064
1065!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1066     !
1067     ! on remet le calendrier a zero
1068     !
1069     IF (raz_date .eq. 1) THEN
1070        itau_phy = 0
1071     ENDIF
1072
1073     !IM cf. AM 081204 BEG
1074     PRINT*,'cycle_diurne3 =',cycle_diurne
1075     !IM cf. AM 081204 END
1076     !
1077     CALL printflag( tabcntr0,radpas,ok_journe, &
1078          ok_instan, ok_region )
1079     !
1080     IF (ABS(dtime-pdtphys).GT.0.001) THEN
1081        WRITE(lunout,*) 'Pas physique n est pas correct',dtime, &
1082             pdtphys
1083        abort_message='Pas physique n est pas correct '
1084        !           call abort_gcm(modname,abort_message,1)
1085        dtime=pdtphys
1086     ENDIF
1087     IF (nlon .NE. klon) THEN
1088        WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,  &
1089             klon
1090        abort_message='nlon et klon ne sont pas coherents'
1091        call abort_gcm(modname,abort_message,1)
1092     ENDIF
1093     IF (nlev .NE. klev) THEN
1094        WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1095             klev
1096        abort_message='nlev et klev ne sont pas coherents'
1097        call abort_gcm(modname,abort_message,1)
1098     ENDIF
1099     !
1100     IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
1101        WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1102        WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1103        abort_message='Nbre d appels au rayonnement insuffisant'
1104        call abort_gcm(modname,abort_message,1)
1105     ENDIF
1106     WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1107     WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
1108          ok_cvl
1109     !
1110     !KE43
1111     ! Initialisation pour la convection de K.E. (sb):
1112     IF (iflag_con.GE.3) THEN
1113
1114        WRITE(lunout,*)"*** Convection de Kerry Emanuel 4.3  "
1115        WRITE(lunout,*) &
1116             "On va utiliser le melange convectif des traceurs qui"
1117        WRITE(lunout,*)"est calcule dans convect4.3"
1118        WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1119
1120        DO i = 1, klon
1121           ema_cbmf(i) = 0.
1122           ema_pcb(i)  = 0.
1123           ema_pct(i)  = 0.
1124           !          ema_workcbmf(i) = 0.
1125        ENDDO
1126        !IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1127        DO i = 1, klon
1128           ibas_con(i) = 1
1129           itop_con(i) = 1
1130        ENDDO
1131        !IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1132        !===============================================================================
1133        !CR:04.12.07: initialisations poches froides
1134        ! Controle de ALE et ALP pour la fermeture convective (jyg)
1135        if (iflag_wake>=1) then
1136           CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1137                ,alp_bl_prescr, ale_bl_prescr)
1138           ! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1139           !        print*,'apres ini_wake iflag_cldth=', iflag_cldth
1140        endif
1141
1142!        do i = 1,klon
1143!           Ale_bl(i)=0.
1144!           Alp_bl(i)=0.
1145!        enddo
1146
1147        !================================================================================
1148        !IM stations CFMIP
1149        nCFMIP=npCFMIP
1150        OPEN(98,file='npCFMIP_param.data',status='old', &
1151             form='formatted',iostat=iostat)
1152        if (iostat == 0) then
1153           READ(98,*,end=998) nCFMIP
1154998        CONTINUE
1155           CLOSE(98)
1156           CONTINUE
1157           IF(nCFMIP.GT.npCFMIP) THEN
1158              print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1159              call abort_gcm("physiq", "", 1)
1160           else
1161              print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1162           ENDIF
1163
1164           !
1165           ALLOCATE(tabCFMIP(nCFMIP))
1166           ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1167           ALLOCATE(tabijGCM(nCFMIP))
1168           ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1169           ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1170           !
1171           ! lecture des nCFMIP stations CFMIP, de leur numero
1172           ! et des coordonnees geographiques lonCFMIP, latCFMIP
1173           !
1174           CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP,  &
1175                lonCFMIP, latCFMIP)
1176           !
1177           ! identification des
1178           ! 1) coordonnees lonGCM, latGCM des points CFMIP dans la grille de LMDZ
1179           ! 2) indices points tabijGCM de la grille physique 1d sur klon points
1180           ! 3) indices iGCM, jGCM de la grille physique 2d
1181           !
1182           CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1183                tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1184           !
1185        else
1186           ALLOCATE(tabijGCM(0))
1187           ALLOCATE(lonGCM(0), latGCM(0))
1188           ALLOCATE(iGCM(0), jGCM(0))
1189        end if
1190     else
1191        ALLOCATE(tabijGCM(0))
1192        ALLOCATE(lonGCM(0), latGCM(0))
1193        ALLOCATE(iGCM(0), jGCM(0))
1194     ENDIF
1195
1196     DO i=1,klon
1197        rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1198     ENDDO
1199
1200     !34EK
1201     IF (ok_orodr) THEN
1202
1203!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1204        ! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1205        ! justement quand ok_orodr = false.
1206        ! ce rugoro est utilise par la couche limite et fait double emploi
1207        ! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1208        !           DO i=1,klon
1209        !             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1210        !           ENDDO
1211!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1212        IF (ok_strato) THEN
1213           CALL SUGWD_strato(klon,klev,paprs,pplay)
1214        ELSE
1215           CALL SUGWD(klon,klev,paprs,pplay)
1216        ENDIF
1217
1218        DO i=1,klon
1219           zuthe(i)=0.
1220           zvthe(i)=0.
1221           if(zstd(i).gt.10.)then
1222              zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1223              zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1224           endif
1225        ENDDO
1226     ENDIF
1227     !
1228     !
1229     lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
1230     WRITE(lunout,*)'La frequence de lecture surface est de ',  &
1231          lmt_pas
1232     !
1233     capemaxcels = 't_max(X)'
1234     t2mincels = 't_min(X)'
1235     t2maxcels = 't_max(X)'
1236     tinst = 'inst(X)'
1237     tave = 'ave(X)'
1238     !IM cf. AM 081204 BEG
1239     write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1240     !IM cf. AM 081204 END
1241     !
1242     !=============================================================
1243     !   Initialisation des sorties
1244     !=============================================================
1245
1246#ifdef CPP_IOIPSL
1247
1248     !$OMP MASTER
1249! FH : if ok_sync=.true. , the time axis is written at each time step
1250! in the output files. Only at the end in the opposite case
1251     ok_sync_omp=.false.
1252     CALL getin('ok_sync',ok_sync_omp)
1253     call phys_output_open(rlon,rlat,nCFMIP,tabijGCM, &
1254          iGCM,jGCM,lonGCM,latGCM, &
1255          jjmp1,nlevSTD,clevSTD,rlevSTD, dtime,ok_veget, &
1256          type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
1257          ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,  &
1258          read_climoz, phys_out_filestations, &
1259          new_aod, aerosol_couple, &
1260          flag_aerosol_strat, pdtphys, paprs, pphis,  &
1261          pplay, lmax_th, ptconv, ptconvth, ivap,  &
1262          d_t, qx, d_qx, zmasse, ok_sync_omp)
1263     !$OMP END MASTER
1264     !$OMP BARRIER
1265     ok_sync=ok_sync_omp
1266
1267     freq_outNMC(1) = ecrit_files(7)
1268     freq_outNMC(2) = ecrit_files(8)
1269     freq_outNMC(3) = ecrit_files(9)
1270     WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1271     WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1272     WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1273
1274     include "ini_histday_seri.h"
1275
1276     include "ini_paramLMDZ_phy.h"
1277
1278#endif
1279     ecrit_reg = ecrit_reg * un_jour
1280     ecrit_tra = ecrit_tra * un_jour
1281
1282     !XXXPB Positionner date0 pour initialisation de ORCHIDEE
1283     date0 = jD_ref
1284     WRITE(*,*) 'physiq date0 : ',date0
1285     !
1286     !
1287     !
1288     ! Prescrire l'ozone dans l'atmosphere
1289     !
1290     !
1291     !c         DO i = 1, klon
1292     !c         DO k = 1, klev
1293     !c            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1294     !c         ENDDO
1295     !c         ENDDO
1296     !
1297     IF (type_trac == 'inca') THEN
1298#ifdef INCA
1299        CALL VTe(VTphysiq)
1300        CALL VTb(VTinca)
1301        calday = REAL(days_elapsed) + jH_cur
1302        WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1303
1304        CALL chemini(  &
1305             rg, &
1306             ra, &
1307             airephy, &
1308             rlat, &
1309             rlon, &
1310             presnivs, &
1311             calday, &
1312             klon, &
1313             nqtot, &
1314             pdtphys, &
1315             annee_ref, &
1316             day_ref,  &
1317             day_ini, &
1318             start_time, &
1319             itau_phy, &
1320             io_lon, &
1321             io_lat)
1322
1323        CALL VTe(VTinca)
1324        CALL VTb(VTphysiq)
1325#endif
1326     END IF
1327     !
1328     !
1329!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1330     ! Nouvelle initialisation pour le rayonnement RRTM
1331!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1332
1333     call iniradia(klon,klev,paprs(1,1:klev+1))
1334
1335     !$omp single
1336     if (read_climoz >= 1) then
1337        call open_climoz(ncid_climoz, press_climoz)
1338     END IF
1339     !$omp end single
1340     !
1341     !IM betaCRF
1342     pfree=70000. !Pa
1343     beta_pbl=1.
1344     beta_free=1.
1345     lon1_beta=-180.
1346     lon2_beta=+180.
1347     lat1_beta=90.
1348     lat2_beta=-90.
1349     mskocean_beta=.FALSE.
1350
1351     OPEN(99,file='beta_crf.data',status='old', &
1352          form='formatted',err=9999)
1353     READ(99,*,end=9998) pfree
1354     READ(99,*,end=9998) beta_pbl
1355     READ(99,*,end=9998) beta_free
1356     READ(99,*,end=9998) lon1_beta
1357     READ(99,*,end=9998) lon2_beta
1358     READ(99,*,end=9998) lat1_beta
1359     READ(99,*,end=9998) lat2_beta
1360     READ(99,*,end=9998) mskocean_beta
13619998 Continue
1362     CLOSE(99)
13639999 Continue
1364     WRITE(*,*)'pfree=',pfree
1365     WRITE(*,*)'beta_pbl=',beta_pbl
1366     WRITE(*,*)'beta_free=',beta_free
1367     WRITE(*,*)'lon1_beta=',lon1_beta
1368     WRITE(*,*)'lon2_beta=',lon2_beta
1369     WRITE(*,*)'lat1_beta=',lat1_beta
1370     WRITE(*,*)'lat2_beta=',lat2_beta
1371     WRITE(*,*)'mskocean_beta=',mskocean_beta
1372  ENDIF
1373  !
1374  !   ****************     Fin  de   IF ( debut  )   ***************
1375  !
1376  !
1377  ! Incrementer le compteur de la physique
1378  !
1379  itap   = itap + 1
1380  !
1381  !
1382  ! Update fraction of the sub-surfaces (pctsrf) and
1383  ! initialize, where a new fraction has appeared, all variables depending
1384  ! on the surface fraction.
1385  !
1386  CALL change_srf_frac(itap, dtime, days_elapsed+1,  &
1387       pctsrf, falb1, falb2, ftsol, ustar, u10m, v10m, pbl_tke)
1388
1389
1390  ! Update time and other variables in Reprobus
1391  IF (type_trac == 'repr') THEN
1392#ifdef REPROBUS
1393     CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1394     print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1395     CALL Rtime(debut)
1396#endif
1397  END IF
1398
1399
1400  ! Tendances bidons pour les processus qui n'affectent pas certaines
1401  ! variables.
1402  du0(:,:)=0.
1403  dv0(:,:)=0.
1404  dt0 = 0.
1405  dq0(:,:)=0.
1406  dql0(:,:)=0.
1407  dqi0(:,:)=0.
1408  !
1409  ! Mettre a zero des variables de sortie (pour securite)
1410  !
1411  DO i = 1, klon
1412     d_ps(i) = 0.0
1413  ENDDO
1414  DO k = 1, klev
1415     DO i = 1, klon
1416        d_t(i,k) = 0.0
1417        d_u(i,k) = 0.0
1418        d_v(i,k) = 0.0
1419     ENDDO
1420  ENDDO
1421  DO iq = 1, nqtot
1422     DO k = 1, klev
1423        DO i = 1, klon
1424           d_qx(i,k,iq) = 0.0
1425        ENDDO
1426     ENDDO
1427  ENDDO
1428  da(:,:)=0.
1429  mp(:,:)=0.
1430  phi(:,:,:)=0.
1431  ! RomP >>>
1432  phi2(:,:,:)=0.
1433  beta_prec_fisrt(:,:)=0.
1434  beta_prec(:,:)=0.
1435  epmlmMm(:,:,:)=0.
1436  eplaMm(:,:)=0.
1437  d1a(:,:)=0.
1438  dam(:,:)=0.
1439  pmflxr=0.
1440  pmflxs=0.
1441  ! RomP <<<
1442
1443  !
1444  ! Ne pas affecter les valeurs entrees de u, v, h, et q
1445  !
1446  DO k = 1, klev
1447     DO i = 1, klon
1448        t_seri(i,k)  = t(i,k)
1449        u_seri(i,k)  = u(i,k)
1450        v_seri(i,k)  = v(i,k)
1451        q_seri(i,k)  = qx(i,k,ivap)
1452        ql_seri(i,k) = qx(i,k,iliq)
1453!CR: ATTENTION, on rajoute la variable glace
1454        if (nqo.eq.2) then
1455           qs_seri(i,k) = 0.
1456        else if (nqo.eq.3) then
1457           qs_seri(i,k) = qx(i,k,isol)
1458        endif
1459     ENDDO
1460  ENDDO
1461  tke0(:,:)=pbl_tke(:,:,is_ave)
1462!CR:Nombre de traceurs de l'eau: nqo
1463!  IF (nqtot.GE.3) THEN
1464   IF (nqtot.GE.(nqo+1)) THEN
1465!     DO iq = 3, nqtot       
1466     DO iq = nqo+1, nqtot 
1467        DO  k = 1, klev
1468           DO  i = 1, klon
1469!              tr_seri(i,k,iq-2) = qx(i,k,iq)
1470              tr_seri(i,k,iq-nqo) = qx(i,k,iq)
1471           ENDDO
1472        ENDDO
1473     ENDDO
1474  ELSE
1475     DO k = 1, klev
1476        DO i = 1, klon
1477           tr_seri(i,k,1) = 0.0
1478        ENDDO
1479     ENDDO
1480  ENDIF
1481  !
1482  DO i = 1, klon
1483     ztsol(i) = 0.
1484  ENDDO
1485  DO nsrf = 1, nbsrf
1486     DO i = 1, klon
1487        ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1488     ENDDO
1489  ENDDO
1490  !IM
1491  IF (ip_ebil_phy.ge.1) THEN
1492     ztit='after dynamic'
1493     CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
1494          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1495          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1496     !     Comme les tendances de la physique sont ajoute dans la dynamique,
1497     !     on devrait avoir que la variation d'entalpie par la dynamique
1498     !     est egale a la variation de la physique au pas de temps precedent.
1499     !     Donc la somme de ces 2 variations devrait etre nulle.
1500     call diagphy(airephy,ztit,ip_ebil_phy &
1501          , zero_v, zero_v, zero_v, zero_v, zero_v &
1502          , zero_v, zero_v, zero_v, ztsol &
1503          , d_h_vcol+d_h_vcol_phy, d_qt, 0. &
1504          , fs_bound, fq_bound )
1505  END IF
1506
1507  ! Diagnostiquer la tendance dynamique
1508  !
1509  IF (ancien_ok) THEN
1510     DO k = 1, klev
1511        DO i = 1, klon
1512           d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1513           d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1514           d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1515           d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1516        ENDDO
1517     ENDDO
1518!!! RomP >>>   td dyn traceur
1519     IF (nqtot.GE.3) THEN
1520        DO iq = 3, nqtot
1521           DO k = 1, klev
1522              DO i = 1, klon
1523                 d_tr_dyn(i,k,iq-2)= &
1524                      (tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime
1525                 !         iiq=niadv(iq)
1526                 !         print*,i,k," d_tr_dyn",d_tr_dyn(i,k,iq-2),"tra:",iq,tname(iiq)
1527              ENDDO
1528           ENDDO
1529        ENDDO
1530     ENDIF
1531!!! RomP <<<
1532  ELSE
1533     DO k = 1, klev
1534        DO i = 1, klon
1535           d_u_dyn(i,k) = 0.0
1536           d_v_dyn(i,k) = 0.0
1537           d_t_dyn(i,k) = 0.0
1538           d_q_dyn(i,k) = 0.0
1539        ENDDO
1540     ENDDO
1541!!! RomP >>>   td dyn traceur
1542     IF (nqtot.GE.3) THEN
1543        DO iq = 3, nqtot
1544           DO k = 1, klev
1545              DO i = 1, klon
1546                 d_tr_dyn(i,k,iq-2)= 0.0
1547              ENDDO
1548           ENDDO
1549        ENDDO
1550     ENDIF
1551!!! RomP <<<
1552     ancien_ok = .TRUE.
1553  ENDIF
1554  !
1555  ! Ajouter le geopotentiel du sol:
1556  !
1557  DO k = 1, klev
1558     DO i = 1, klon
1559        zphi(i,k) = pphi(i,k) + pphis(i)
1560     ENDDO
1561  ENDDO
1562  !
1563  ! Verifier les temperatures
1564  !
1565  !IM BEG
1566  IF (check) THEN
1567     amn=MIN(ftsol(1,is_ter),1000.)
1568     amx=MAX(ftsol(1,is_ter),-1000.)
1569     DO i=2, klon
1570        amn=MIN(ftsol(i,is_ter),amn)
1571        amx=MAX(ftsol(i,is_ter),amx)
1572     ENDDO
1573     !
1574     PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1575  ENDIF !(check) THEN
1576  !IM END
1577  !
1578  CALL hgardfou(t_seri,ftsol,'debutphy')
1579  !
1580  !IM BEG
1581  IF (check) THEN
1582     amn=MIN(ftsol(1,is_ter),1000.)
1583     amx=MAX(ftsol(1,is_ter),-1000.)
1584     DO i=2, klon
1585        amn=MIN(ftsol(i,is_ter),amn)
1586        amx=MAX(ftsol(i,is_ter),amx)
1587     ENDDO
1588     !
1589     PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1590  ENDIF !(check) THEN
1591  !IM END
1592  !
1593  ! Mettre en action les conditions aux limites (albedo, sst, etc.).
1594  ! Prescrire l'ozone et calculer l'albedo sur l'ocean.
1595  !
1596  if (read_climoz >= 1) then
1597     ! Ozone from a file
1598     ! Update required ozone index:
1599     ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1
1600     if (ro3i == 361) ro3i = 360
1601     ! (This should never occur, except perhaps because of roundup
1602     ! error. See documentation.)
1603     if (ro3i /= co3i) then
1604        ! Update ozone field:
1605        if (read_climoz == 1) then
1606           call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
1607                press_in_edg=press_climoz, paprs=paprs, v3=wo)
1608        else
1609           ! read_climoz == 2
1610           call regr_pr_av(ncid_climoz, (/"tro3         ", "tro3_daylight"/), &
1611                julien=ro3i, press_in_edg=press_climoz, paprs=paprs, v3=wo)
1612        end if
1613        ! Convert from mole fraction of ozone to column density of ozone in a
1614        ! cell, in kDU:
1615        forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
1616             * zmasse / dobson_u / 1e3
1617        ! (By regridding ozone values for LMDZ only once every 360th of
1618        ! year, we have already neglected the variation of pressure in one
1619        ! 360th of year. So do not recompute "wo" at each time step even if
1620        ! "zmasse" changes a little.)
1621        co3i = ro3i
1622     end if
1623  ELSEIF (MOD(itap-1,lmt_pas) == 0) THEN
1624     ! Once per day, update ozone from Royer:
1625
1626     IF (solarlong0<-999.) then
1627        ! Generic case with evolvoing season
1628        zzz=real(days_elapsed+1)
1629     ELSE IF (abs(solarlong0-1000.)<1.e-4) then
1630        ! Particular case with annual mean insolation
1631        zzz=real(90) ! could be revisited
1632        IF (read_climoz/=-1) THEN
1633           abort_message ='read_climoz=-1 is recommended when solarlong0=1000.'
1634           CALL abort_gcm (modname,abort_message,1)
1635        ENDIF
1636     ELSE
1637        ! Case where the season is imposed with solarlong0
1638        zzz=real(90) ! could be revisited
1639     ENDIF
1640     wo(:,:,1)=ozonecm(rlat, paprs,read_climoz,rjour=zzz)
1641  ENDIF
1642  !
1643  ! Re-evaporer l'eau liquide nuageuse
1644  !
1645  DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
1646     DO i = 1, klon
1647        zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1648        !jyg<
1649        !      Attention : Arnaud a propose des formules completement differentes
1650        !                  A verifier !!!
1651        zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1652        IF (iflag_ice_thermo .EQ. 0) THEN
1653           zlsdcp=zlvdcp
1654        ENDIF
1655        !>jyg
1656     
1657        if (iflag_ice_thermo.eq.0) then   
1658!pas necessaire a priori
1659
1660        zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1661        zb = MAX(0.0,ql_seri(i,k))
1662        za = - MAX(0.0,ql_seri(i,k)) &
1663             * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
1664        t_seri(i,k) = t_seri(i,k) + za
1665        q_seri(i,k) = q_seri(i,k) + zb
1666        ql_seri(i,k) = 0.0
1667        d_t_eva(i,k) = za
1668        d_q_eva(i,k) = zb
1669
1670        else
1671
1672!CR: on ré-évapore eau liquide et glace
1673
1674!        zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1675!        zb = MAX(0.0,ql_seri(i,k))
1676!        za = - MAX(0.0,ql_seri(i,k)) &
1677!             * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
1678        zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
1679        za = - MAX(0.0,ql_seri(i,k))*zlvdcp &
1680             - MAX(0.0,qs_seri(i,k))*zlsdcp
1681        t_seri(i,k) = t_seri(i,k) + za
1682        q_seri(i,k) = q_seri(i,k) + zb
1683        ql_seri(i,k) = 0.0
1684!on évapore la glace
1685        qs_seri(i,k) = 0.0
1686        d_t_eva(i,k) = za
1687        d_q_eva(i,k) = zb
1688        endif
1689
1690     ENDDO
1691  ENDDO
1692  !IM
1693  IF (ip_ebil_phy.ge.2) THEN
1694     ztit='after reevap'
1695     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime &
1696          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1697          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1698     call diagphy(airephy,ztit,ip_ebil_phy &
1699          , zero_v, zero_v, zero_v, zero_v, zero_v &
1700          , zero_v, zero_v, zero_v, ztsol &
1701          , d_h_vcol, d_qt, d_ec &
1702          , fs_bound, fq_bound )
1703     !
1704  END IF
1705
1706  !
1707  !=========================================================================
1708  ! Calculs de l'orbite.
1709  ! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1710  ! doit donc etre plac\'e avant radlwsw et pbl_surface
1711
1712!!!   jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1713  call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1714  day_since_equinox = (jD_cur + jH_cur) - jD_eq
1715  !
1716  !   choix entre calcul de la longitude solaire vraie ou valeur fixee a
1717  !   solarlong0
1718  if (solarlong0<-999.) then
1719     if (new_orbit) then
1720        ! calcul selon la routine utilisee pour les planetes
1721        call solarlong(day_since_equinox, zlongi, dist)
1722     else
1723        ! calcul selon la routine utilisee pour l'AR4
1724        CALL orbite(REAL(days_elapsed+1),zlongi,dist)
1725     endif
1726  else
1727     zlongi=solarlong0  ! longitude solaire vraie
1728     dist=1.            ! distance au soleil / moyenne
1729  endif
1730  if(prt_level.ge.1)                                                &
1731       write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1732
1733
1734!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1735  ! Calcul de l'ensoleillement :
1736  ! ============================
1737  ! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
1738  ! l'annee a partir d'une formule analytique.
1739  ! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
1740  ! non nul aux poles.
1741  IF (abs(solarlong0-1000.)<1.e-4) then
1742     call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
1743  ELSE
1744     !  Avec ou sans cycle diurne
1745     IF (cycle_diurne) THEN
1746        zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1747        CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
1748     ELSE
1749        CALL angle(zlongi, rlat, fract, rmu0)
1750     ENDIF
1751  ENDIF
1752
1753  ! AI Janv 2014
1754  do i = 1, klon
1755     if (fract(i).le.0.) then
1756        JrNt(i)=0.
1757     else
1758        JrNt(i)=1.
1759     endif
1760  enddo
1761
1762  if (mydebug) then
1763     call writefield_phy('u_seri',u_seri,llm)
1764     call writefield_phy('v_seri',v_seri,llm)
1765     call writefield_phy('t_seri',t_seri,llm)
1766     call writefield_phy('q_seri',q_seri,llm)
1767  endif
1768
1769  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1770  ! Appel au pbl_surface : Planetary Boudary Layer et Surface
1771  ! Cela implique tous les interactions des sous-surfaces et la partie diffusion
1772  ! turbulent du couche limit.
1773  !
1774  ! Certains varibales de sorties de pbl_surface sont utiliser que pour
1775  ! ecriture des fihiers hist_XXXX.nc, ces sont :
1776  !   qsol,      zq2m,      s_pblh,  s_lcl,
1777  !   s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
1778  !   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
1779  !   zxrugs,    zu10m,     zv10m,   fder,
1780  !   zxqsurf,   rh2m,      zxfluxu, zxfluxv,
1781  !   frugs,     agesno,    fsollw,  fsolsw,
1782  !   d_ts,      fevap,     fluxlat, t2m,
1783  !   wfbils,    wfbilo,    fluxt,   fluxu, fluxv,
1784  !
1785  ! Certains ne sont pas utiliser du tout :
1786  !   dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1787  !
1788
1789  ! Calcul de l'humidite de saturation au niveau du sol
1790
1791
1792
1793  if (iflag_pbl/=0) then
1794
1795!jyg+nrlmd<
1796      IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
1797        print *,'debut du splitting de la PBL'
1798      ENDIF
1799!!!
1800!=================================================================
1801!         PROVISOIRE : DECOUPLAGE PBL/WAKE
1802!         --------------------------------
1803!
1804!!      wake_deltat_sav(:,:)=wake_deltat(:,:)
1805!!      wake_deltaq_sav(:,:)=wake_deltaq(:,:)
1806!!      wake_deltat(:,:)=0.
1807!!      wake_deltaq(:,:)=0.
1808!=================================================================
1809!>jyg+nrlmd
1810!
1811     CALL pbl_surface(  &
1812          dtime,     date0,     itap,    days_elapsed+1, &
1813          debut,     lafin, &
1814          rlon,      rlat,      rugoro,  rmu0,      &
1815          zsig,      sollwdown, pphi,    cldt,      &
1816          rain_fall, snow_fall, solsw,   sollw,     &
1817          t_seri,    q_seri,    u_seri,  v_seri,    &
1818!nrlmd+jyg<
1819          wake_deltat, wake_deltaq, wake_cstar, wake_s, &
1820!>nrlmd+jyg
1821          pplay,     paprs,     pctsrf,             &
1822          ftsol,falb1,falb2,ustar,u10m,v10m,wstar,  &
1823          cdragh,    cdragm,  u1,    v1,            &
1824          albsol1,   albsol2,   sens,    evap,      &
1825          albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
1826          zxtsol,    zxfluxlat, zt2m,    qsat2m,  &
1827          d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss, &
1828!nrlmd<
1829  !jyg<
1830          d_t_vdf_w, d_q_vdf_w, &
1831          d_t_vdf_x, d_q_vdf_x, &
1832          sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
1833  !>jyg
1834          delta_tsurf,wake_dens, &
1835          cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
1836          kh,kh_x,kh_w, &
1837!>nrlmd
1838          coefh(1:klon,1:klev,1:nbsrf+1),     coefm(1:klon,1:klev,1:nbsrf+1), &
1839          slab_wfbils,                 &
1840          qsol,      zq2m,      s_pblh,  s_lcl, &
1841!jyg<
1842          s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
1843!>jyg
1844          s_capCL,   s_oliqCL,  s_cteiCL,s_pblT, &
1845          s_therm,   s_trmb1,   s_trmb2, s_trmb3, &
1846          zxrugs,    zustar, zu10m,     zv10m,   fder, &
1847          zxqsurf,   rh2m,      zxfluxu, zxfluxv, &
1848          frugs,     agesno,    fsollw,  fsolsw, &
1849          d_ts,      fevap,     fluxlat, t2m, &
1850          wfbils,    wfbilo,    fluxt,   fluxu,  fluxv, &
1851          dsens,     devap,     zxsnow, &
1852          zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke, &
1853!nrlmd+jyg<
1854          wake_delta_pbl_TKE &
1855!>nrlmd+jyg
1856                      )
1857!
1858!=================================================================
1859!         PROVISOIRE : DECOUPLAGE PBL/WAKE
1860!         --------------------------------
1861!
1862!!      wake_deltat(:,:)=wake_deltat_sav(:,:)
1863!!      wake_deltaq(:,:)=wake_deltaq_sav(:,:)
1864!=================================================================
1865!
1866!  Add turbulent diffusion tendency to the wake difference variables
1867    IF (mod(iflag_pbl_split,2) .NE. 0) THEN
1868     wake_deltat(:,:) = wake_deltat(:,:) + (d_t_vdf_w(:,:)-d_t_vdf_x(:,:))
1869     wake_deltaq(:,:) = wake_deltaq(:,:) + (d_q_vdf_w(:,:)-d_q_vdf_x(:,:))
1870    ENDIF
1871
1872
1873     !---------------------------------------------------------------------
1874     ! ajout des tendances de la diffusion turbulente
1875     IF (klon_glo==1) THEN
1876        CALL add_pbl_tend &
1877             (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf')
1878     ELSE
1879        CALL add_phys_tend &
1880             (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf')
1881     ENDIF
1882     !--------------------------------------------------------------------
1883
1884     if (mydebug) then
1885        call writefield_phy('u_seri',u_seri,llm)
1886        call writefield_phy('v_seri',v_seri,llm)
1887        call writefield_phy('t_seri',t_seri,llm)
1888        call writefield_phy('q_seri',q_seri,llm)
1889     endif
1890
1891     CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
1892          t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
1893
1894
1895     IF (ip_ebil_phy.ge.2) THEN
1896        ztit='after surface_main'
1897        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
1898             , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1899             , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1900        call diagphy(airephy,ztit,ip_ebil_phy &
1901             , zero_v, zero_v, zero_v, zero_v, sens &
1902             , evap  , zero_v, zero_v, ztsol &
1903             , d_h_vcol, d_qt, d_ec &
1904             , fs_bound, fq_bound )
1905     END IF
1906
1907  ENDIF
1908  ! =================================================================== c
1909  !   Calcul de Qsat
1910
1911  DO k = 1, klev
1912     DO i = 1, klon
1913        zx_t = t_seri(i,k)
1914        IF (thermcep) THEN
1915           zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1916           zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1917           zx_qs  = MIN(0.5,zx_qs)
1918           zcor   = 1./(1.-retv*zx_qs)
1919           zx_qs  = zx_qs*zcor
1920        ELSE
1921           IF (zx_t.LT.t_coup) THEN
1922              zx_qs = qsats(zx_t)/pplay(i,k)
1923           ELSE
1924              zx_qs = qsatl(zx_t)/pplay(i,k)
1925           ENDIF
1926        ENDIF
1927        zqsat(i,k)=zx_qs
1928     ENDDO
1929  ENDDO
1930
1931  if (prt_level.ge.1) then
1932     write(lunout,*) 'L   qsat (g/kg) avant clouds_gno'
1933     write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1934  endif
1935  !
1936  ! Appeler la convection (au choix)
1937  !
1938  DO k = 1, klev
1939     DO i = 1, klon
1940        conv_q(i,k) = d_q_dyn(i,k)  &
1941             + d_q_vdf(i,k)/dtime
1942        conv_t(i,k) = d_t_dyn(i,k)  &
1943             + d_t_vdf(i,k)/dtime
1944     ENDDO
1945  ENDDO
1946  IF (check) THEN
1947     za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1948     WRITE(lunout,*) "avantcon=", za
1949  ENDIF
1950  zx_ajustq = .FALSE.
1951  IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1952  IF (zx_ajustq) THEN
1953     DO i = 1, klon
1954        z_avant(i) = 0.0
1955     ENDDO
1956     DO k = 1, klev
1957        DO i = 1, klon
1958           z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
1959                *(paprs(i,k)-paprs(i,k+1))/RG
1960        ENDDO
1961     ENDDO
1962  ENDIF
1963
1964  ! Calcule de vitesse verticale a partir de flux de masse verticale
1965  DO k = 1, klev
1966     DO i = 1, klon
1967        omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
1968     END DO
1969  END DO
1970  if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
1971       omega(igout, :)
1972
1973  IF (iflag_con.EQ.1) THEN
1974     abort_message ='reactiver le call conlmd dans physiq.F'
1975     CALL abort_gcm (modname,abort_message,1)
1976     !     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1977     !    .             d_t_con, d_q_con,
1978     !    .             rain_con, snow_con, ibas_con, itop_con)
1979  ELSE IF (iflag_con.EQ.2) THEN
1980     CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
1981          conv_t, conv_q, -evap, omega, &
1982          d_t_con, d_q_con, rain_con, snow_con, &
1983          pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
1984          kcbot, kctop, kdtop, pmflxr, pmflxs)
1985     d_u_con = 0.
1986     d_v_con = 0.
1987
1988     WHERE (rain_con < 0.) rain_con = 0.
1989     WHERE (snow_con < 0.) snow_con = 0.
1990     DO i = 1, klon
1991        ibas_con(i) = klev+1 - kcbot(i)
1992        itop_con(i) = klev+1 - kctop(i)
1993     ENDDO
1994  ELSE IF (iflag_con.GE.3) THEN
1995     ! nb of tracers for the KE convection:
1996     ! MAF la partie traceurs est faite dans phytrac
1997     ! on met ntra=1 pour limiter les appels mais on peut
1998     ! supprimer les calculs / ftra.
1999     ntra = 1
2000
2001     !=========================================================================
2002     !ajout pour la parametrisation des poches froides: calcul de
2003     !t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2004     do k=1,klev
2005        do i=1,klon
2006           if (iflag_wake>=1) then
2007              t_wake(i,k) = t_seri(i,k) &
2008                   +(1-wake_s(i))*wake_deltat(i,k)
2009              q_wake(i,k) = q_seri(i,k) &
2010                   +(1-wake_s(i))*wake_deltaq(i,k)
2011              t_undi(i,k) = t_seri(i,k) &
2012                   -wake_s(i)*wake_deltat(i,k)
2013              q_undi(i,k) = q_seri(i,k) &
2014                   -wake_s(i)*wake_deltaq(i,k)
2015           else
2016              t_wake(i,k) = t_seri(i,k)
2017              q_wake(i,k) = q_seri(i,k)
2018              t_undi(i,k) = t_seri(i,k)
2019              q_undi(i,k) = q_seri(i,k)
2020           endif
2021        enddo
2022     enddo
2023
2024     ! Calcul de l'energie disponible ALE (J/kg) et de la puissance
2025     ! disponible ALP (W/m2) pour le soulevement des particules dans
2026     ! le modele convectif
2027     !
2028     do i = 1,klon
2029        ALE(i) = 0.
2030        ALP(i) = 0.
2031     enddo
2032     !
2033     !calcul de ale_wake et alp_wake
2034     if (iflag_wake>=1) then
2035        if (itap .le. it_wape_prescr) then
2036           do i = 1,klon
2037              ale_wake(i) = wape_prescr
2038              alp_wake(i) = fip_prescr
2039           enddo
2040        else
2041           do i = 1,klon
2042              !jyg  ALE=WAPE au lieu de ALE = 1/2 Cstar**2
2043              !cc           ale_wake(i) = 0.5*wake_cstar(i)**2
2044              ale_wake(i) = wake_pe(i)
2045              alp_wake(i) = wake_fip(i)
2046           enddo
2047        endif
2048     else
2049        do i = 1,klon
2050           ale_wake(i) = 0.
2051           alp_wake(i) = 0.
2052        enddo
2053     endif
2054     !combinaison avec ale et alp de couche limite: constantes si pas
2055     !de couplage, valeurs calculees dans le thermique sinon
2056     if (iflag_coupl.eq.0) then
2057        if (debut.and.prt_level.gt.9) &
2058             WRITE(lunout,*)'ALE et ALP imposes'
2059        do i = 1,klon
2060           !on ne couple que ale
2061           !           ALE(i) = max(ale_wake(i),Ale_bl(i))
2062           ALE(i) = max(ale_wake(i),ale_bl_prescr)
2063           !on ne couple que alp
2064           !           ALP(i) = alp_wake(i) + Alp_bl(i)
2065           ALP(i) = alp_wake(i) + alp_bl_prescr
2066        enddo
2067     else
2068        IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2069        !         do i = 1,klon
2070        !             ALE(i) = max(ale_wake(i),Ale_bl(i))
2071        ! avant        ALP(i) = alp_wake(i) + Alp_bl(i)
2072        !             ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2073        !         write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2074        !         write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2075        !         enddo
2076
2077        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2078        ! Modif FH 2010/04/27. Sans doute temporaire.
2079        ! Deux options pour le alp_offset : constant si >?? 0 ou
2080        ! proportionnel ??a w si <0
2081        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2082! Estimation d'une vitesse verticale effective pour ALP
2083        www(1:klon)=0.
2084        do k=2,klev-1
2085           do i=1,klon
2086              www(i)=max(www(i),-omega(i,k)*RD*t_seri(i,k)/(RG*paprs(i,k)) &
2087&                    *zw2(i,k)*zw2(i,k))
2088!             if (paprs(i,k)>pbase(i)) then
2089! calcul approche de la vitesse verticale en m/s
2090!                www(i)=max(www(i),-omega(i,k)*RD*temp(i,k)/(RG*paprs(i,k))
2091!             endif
2092!   Le 0.1 est en gros H / ps = 1e5 / 1e4
2093           enddo
2094        enddo
2095        do i=1,klon
2096           if (www(i)>0. .and. ale_bl(i)>0. ) www(i)=www(i)/ale_bl(i)
2097        enddo
2098
2099
2100        do i = 1,klon
2101           ALE(i) = max(ale_wake(i),Ale_bl(i))
2102           !cc nrlmd le 10/04/2012----------Stochastic triggering--------------
2103           if (iflag_trig_bl.ge.1) then
2104              ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
2105           endif
2106           !cc fin nrlmd le 10/04/2012
2107           if (alp_offset>=0.) then
2108              ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2109           else
2110
2111!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2112!                                _                  _
2113! Ajout d'une composante 3 * A * w w'2 a  w'3  avec w=www : w max sous pbase
2114!       ou A est la fraction couverte par les ascendances w'
2115!       on utilise le fait que A * w'3 = ALP
2116!       et donc A * w'2 ~ ALP / sqrt(ALE)  (on ajoute 0.1 pour les
2117!       singularites)
2118             ALP(i)=alp_wake(i)*(1.+3.*www(i)/( sqrt(ale_wake(i))+0.1) ) &
2119  &                +alp_bl(i)  *(1.+3.*www(i)/( sqrt(ale_bl(i))  +0.1) )
2120!             ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2121!             if (alp(i)<0.) then
2122!                print*,'ALP ',alp(i),alp_wake(i) &
2123!                     ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2124!             endif
2125           endif
2126        enddo
2127!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2128
2129     endif
2130     do i=1,klon
2131        if (alp(i)>alp_max) then
2132           IF(prt_level>9)WRITE(lunout,*)                             &
2133                'WARNING SUPER ALP (seuil=',alp_max, &
2134                '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2135           alp(i)=alp_max
2136        endif
2137        if (ale(i)>ale_max) then
2138           IF(prt_level>9)WRITE(lunout,*)                             &
2139                'WARNING SUPER ALE (seuil=',ale_max, &
2140                '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2141           ale(i)=ale_max
2142        endif
2143     enddo
2144
2145     !fin calcul ale et alp
2146     !=======================================================================
2147
2148
2149     ! sb, oct02:
2150     ! Schema de convection modularise et vectorise:
2151     ! (driver commun aux versions 3 et 4)
2152     !
2153     IF (ok_cvl) THEN ! new driver for convectL
2154
2155        IF (type_trac == 'repr') THEN
2156           nbtr_tmp=ntra
2157        ELSE
2158           nbtr_tmp=nbtr
2159        END IF
2160        !jyg   iflag_con est dans clesphys
2161        !c          CALL concvl (iflag_con,iflag_clos,
2162        CALL concvl (iflag_clos, &
2163             dtime,paprs,pplay,t_undi,q_undi, &
2164             t_wake,q_wake,wake_s, &
2165             u_seri,v_seri,tr_seri,nbtr_tmp, &
2166             ALE,ALP, &
2167             sig1,w01, &
2168             d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2169             rain_con, snow_con, ibas_con, itop_con, sigd, &
2170             ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0, &
2171             Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2172             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2173             ! RomP >>>
2174             !!     .        pmflxr,pmflxs,da,phi,mp,
2175             !!     .        ftd,fqd,lalim_conv,wght_th)
2176             pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
2177             ftd,fqd,lalim_conv,wght_th, &
2178             ev, ep,epmlmMm,eplaMm, &
2179             wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
2180             tau_cld_cv,coefw_cld_cv)
2181        ! RomP <<<
2182
2183        !IM begin
2184        !       print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2185        !    .dnwd0(1,1),ftd(1,1),fqd(1,1)
2186        !IM end
2187        !IM cf. FH
2188        clwcon0=qcondc
2189        pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2190
2191        do i = 1, klon
2192           if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2193        enddo
2194
2195     ELSE ! ok_cvl
2196
2197        ! MAF conema3 ne contient pas les traceurs
2198        CALL conema3 (dtime, &
2199             paprs,pplay,t_seri,q_seri, &
2200             u_seri,v_seri,tr_seri,ntra, &
2201             sig1,w01, &
2202             d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2203             rain_con, snow_con, ibas_con, itop_con, &
2204             upwd,dnwd,dnwd0,bas,top, &
2205             Ma,cape,tvp,rflag, &
2206             pbase &
2207             ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
2208             ,clwcon0)
2209
2210     ENDIF ! ok_cvl
2211
2212     !
2213     ! Correction precip
2214     rain_con = rain_con * cvl_corr
2215     snow_con = snow_con * cvl_corr
2216     !
2217
2218     IF (.NOT. ok_gust) THEN
2219        do i = 1, klon
2220           wd(i)=0.0
2221        enddo
2222     ENDIF
2223
2224     ! =================================================================== c
2225     ! Calcul des proprietes des nuages convectifs
2226     !
2227
2228     !   calcul des proprietes des nuages convectifs
2229     clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2230     IF (iflag_cld_cv <= 1) THEN
2231     call clouds_gno &
2232          (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2233     ELSE
2234     call clouds_bigauss &
2235          (klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
2236     ENDIF
2237
2238
2239     ! =================================================================== c
2240
2241     DO i = 1, klon
2242        itop_con(i) = min(max(itop_con(i),1),klev)
2243        ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2244     ENDDO
2245
2246     DO i = 1, klon
2247        ema_pcb(i)  = paprs(i,ibas_con(i))
2248     ENDDO
2249     DO i = 1, klon
2250        ! L'idicage de itop_con peut cacher un pb potentiel
2251        ! FH sous la dictee de JYG, CR
2252        ema_pct(i)  = paprs(i,itop_con(i)+1)
2253
2254        if (itop_con(i).gt.klev-3) then
2255           if(prt_level >= 9) then
2256              write(lunout,*)'La convection monte trop haut '
2257              write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2258           endif
2259        endif
2260     ENDDO
2261  ELSE IF (iflag_con.eq.0) THEN
2262     write(lunout,*) 'On n appelle pas la convection'
2263     clwcon0=0.
2264     rnebcon0=0.
2265     d_t_con=0.
2266     d_q_con=0.
2267     d_u_con=0.
2268     d_v_con=0.
2269     rain_con=0.
2270     snow_con=0.
2271     bas=1
2272     top=1
2273  ELSE
2274     WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2275     call abort_gcm("physiq", "", 1)
2276  ENDIF
2277
2278  !     CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2279  !    .              d_u_con, d_v_con)
2280
2281  CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
2282       'convection')
2283  !----------------------------------------------------------------------------
2284
2285  if (mydebug) then
2286     call writefield_phy('u_seri',u_seri,llm)
2287     call writefield_phy('v_seri',v_seri,llm)
2288     call writefield_phy('t_seri',t_seri,llm)
2289     call writefield_phy('q_seri',q_seri,llm)
2290  endif
2291
2292  !IM
2293  IF (ip_ebil_phy.ge.2) THEN
2294     ztit='after convect'
2295     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2296          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2297          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2298     call diagphy(airephy,ztit,ip_ebil_phy &
2299          , zero_v, zero_v, zero_v, zero_v, zero_v &
2300          , zero_v, rain_con, snow_con, ztsol &
2301          , d_h_vcol, d_qt, d_ec &
2302          , fs_bound, fq_bound )
2303  END IF
2304  !
2305  IF (check) THEN
2306     za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2307     WRITE(lunout,*)"aprescon=", za
2308     zx_t = 0.0
2309     za = 0.0
2310     DO i = 1, klon
2311        za = za + airephy(i)/REAL(klon)
2312        zx_t = zx_t + (rain_con(i)+ &
2313             snow_con(i))*airephy(i)/REAL(klon)
2314     ENDDO
2315     zx_t = zx_t/za*dtime
2316     WRITE(lunout,*)"Precip=", zx_t
2317  ENDIF
2318  IF (zx_ajustq) THEN
2319     DO i = 1, klon
2320        z_apres(i) = 0.0
2321     ENDDO
2322     DO k = 1, klev
2323        DO i = 1, klon
2324           z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
2325                *(paprs(i,k)-paprs(i,k+1))/RG
2326        ENDDO
2327     ENDDO
2328     DO i = 1, klon
2329        z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
2330             /z_apres(i)
2331     ENDDO
2332     DO k = 1, klev
2333        DO i = 1, klon
2334           IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
2335                z_factor(i).LT.(1.0-1.0E-08)) THEN
2336              q_seri(i,k) = q_seri(i,k) * z_factor(i)
2337           ENDIF
2338        ENDDO
2339     ENDDO
2340  ENDIF
2341  zx_ajustq=.FALSE.
2342
2343  !
2344  !=============================================================================
2345  !RR:Evolution de la poche froide: on ne fait pas de separation wake/env
2346  !pour la couche limite diffuse pour l instant
2347  !
2348  !
2349  !!! nrlmd le 22/03/2011---Si on met les poches hors des thermiques il faut rajouter cette
2350  !------------------------- tendance calculée hors des poches froides
2351  !
2352  if (iflag_wake>=1) then
2353     DO k=1,klev
2354        DO i=1,klon
2355           dt_dwn(i,k)  = ftd(i,k)
2356           dq_dwn(i,k)  = fqd(i,k)
2357           M_dwn(i,k)   = dnwd0(i,k)
2358           M_up(i,k)    = upwd(i,k)
2359           dt_a(i,k)    = d_t_con(i,k)/dtime - ftd(i,k)
2360           dq_a(i,k)    = d_q_con(i,k)/dtime - fqd(i,k)
2361        ENDDO
2362     ENDDO
2363!nrlmd+jyg<
2364     DO k=1,klev
2365        DO i=1,klon
2366          wdt_PBL(i,k) =  0.
2367          wdq_PBL(i,k) =  0.
2368          udt_PBL(i,k) =  0.
2369          udq_PBL(i,k) =  0.
2370        ENDDO
2371     ENDDO
2372!
2373     IF (mod(iflag_pbl_split,2) .EQ. 1) THEN
2374       DO k=1,klev
2375        DO i=1,klon
2376       wdt_PBL(i,k) = wdt_PBL(i,k) + d_t_vdf_w(i,k)/dtime
2377       wdq_PBL(i,k) = wdq_PBL(i,k) + d_q_vdf_w(i,k)/dtime
2378       udt_PBL(i,k) = udt_PBL(i,k) + d_t_vdf_x(i,k)/dtime
2379       udq_PBL(i,k) = udq_PBL(i,k) + d_q_vdf_x(i,k)/dtime
2380!!        dt_dwn(i,k)  = dt_dwn(i,k) + d_t_vdf_w(i,k)/dtime
2381!!        dq_dwn(i,k)  = dq_dwn(i,k) + d_q_vdf_w(i,k)/dtime
2382!!        dt_a  (i,k)    = dt_a(i,k) + d_t_vdf_x(i,k)/dtime
2383!!        dq_a  (i,k)    = dq_a(i,k) + d_q_vdf_x(i,k)/dtime
2384        ENDDO
2385       ENDDO
2386      ENDIF
2387      IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2388       DO k=1,klev
2389        DO i=1,klon
2390!!        dt_dwn(i,k)  = dt_dwn(i,k) + 0.
2391!!        dq_dwn(i,k)  = dq_dwn(i,k) + 0.
2392!!        dt_a(i,k)   = dt_a(i,k)   + d_t_ajs(i,k)/dtime
2393!!        dq_a(i,k)   = dq_a(i,k)   + d_q_ajs(i,k)/dtime
2394        udt_PBL(i,k)   = udt_PBL(i,k)   + d_t_ajs(i,k)/dtime
2395        udq_PBL(i,k)   = udq_PBL(i,k)   + d_q_ajs(i,k)/dtime
2396        ENDDO
2397       ENDDO
2398      ENDIF
2399!>nrlmd+jyg
2400
2401     IF (iflag_wake==2) THEN
2402        ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2403        DO k = 1,klev
2404           dt_dwn(:,k)= dt_dwn(:,k)+ &
2405                ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
2406           dq_dwn(:,k)= dq_dwn(:,k)+ &
2407                ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
2408        ENDDO
2409     ELSEIF (iflag_wake==3) THEN
2410        ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2411        DO k = 1,klev
2412           DO i=1,klon
2413              IF (rneb(i,k)==0.) THEN
2414! On ne tient compte des tendances qu'en dehors des nuages (c'est �|  dire
2415! a priri dans une region ou l'eau se reevapore).
2416                dt_dwn(i,k)= dt_dwn(i,k)+ &
2417                ok_wk_lsp(i)*d_t_lsc(i,k)/dtime
2418                dq_dwn(i,k)= dq_dwn(i,k)+ &
2419                ok_wk_lsp(i)*d_q_lsc(i,k)/dtime
2420              ENDIF
2421           ENDDO
2422        ENDDO
2423     ENDIF
2424
2425     !
2426     !calcul caracteristiques de la poche froide
2427     call calWAKE (paprs,pplay,dtime &
2428          ,t_seri,q_seri,omega &
2429          ,dt_dwn,dq_dwn,M_dwn,M_up &
2430          ,dt_a,dq_a,sigd &
2431          ,wdt_PBL,wdq_PBL &
2432          ,udt_PBL,udq_PBL &
2433          ,wake_deltat,wake_deltaq,wake_dth &
2434          ,wake_h,wake_s,wake_dens &
2435          ,wake_pe,wake_fip,wake_gfl &
2436          ,dt_wake,dq_wake &
2437          ,wake_k, t_undi,q_undi &
2438          ,wake_omgbdth,wake_dp_omgb &
2439          ,wake_dtKE,wake_dqKE &
2440          ,wake_dtPBL,wake_dqPBL &
2441          ,wake_omg,wake_dp_deltomg &
2442          ,wake_spread,wake_Cstar,wake_d_deltat_gw &
2443          ,wake_ddeltat,wake_ddeltaq)
2444     !
2445     !-------------------------------------------------------------------------
2446     ! ajout des tendances des poches froides
2447     ! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2448     ! coherent avec les autres d_t_...
2449     d_t_wake(:,:)=dt_wake(:,:)*dtime
2450     d_q_wake(:,:)=dq_wake(:,:)*dtime
2451     CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake')
2452     !------------------------------------------------------------------------
2453
2454  endif  ! (iflag_wake>=1)
2455  !
2456  !===================================================================
2457  !JYG
2458  IF (ip_ebil_phy.ge.2) THEN
2459     ztit='after wake'
2460     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2461          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2462          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2463     call diagphy(airephy,ztit,ip_ebil_phy &
2464          , zero_v, zero_v, zero_v, zero_v, zero_v &
2465          , zero_v, zero_v, zero_v, ztsol &
2466          , d_h_vcol, d_qt, d_ec &
2467          , fs_bound, fq_bound )
2468  END IF
2469
2470  !      print*,'apres callwake iflag_cldth=', iflag_cldth
2471  !
2472  !===================================================================
2473  ! Convection seche (thermiques ou ajustement)
2474  !===================================================================
2475  !
2476  call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
2477       ,seuil_inversion,weak_inversion,dthmin)
2478
2479
2480
2481  d_t_ajsb(:,:)=0.
2482  d_q_ajsb(:,:)=0.
2483  d_t_ajs(:,:)=0.
2484  d_u_ajs(:,:)=0.
2485  d_v_ajs(:,:)=0.
2486  d_q_ajs(:,:)=0.
2487  clwcon0th(:,:)=0.
2488  !
2489  !      fm_therm(:,:)=0.
2490  !      entr_therm(:,:)=0.
2491  !      detr_therm(:,:)=0.
2492  !
2493  IF(prt_level>9)WRITE(lunout,*) &
2494       'AVANT LA CONVECTION SECHE , iflag_thermals=' &
2495       ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
2496  if(iflag_thermals<0) then
2497     !  Rien
2498     !  ====
2499     IF(prt_level>9)WRITE(lunout,*)'pas de convection seche'
2500
2501
2502  else
2503
2504     !  Thermiques
2505     !  ==========
2506     IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
2507          ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
2508
2509
2510     !cc nrlmd le 10/04/2012
2511     DO k=1,klev+1
2512        DO i=1,klon
2513           pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
2514           pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
2515           pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
2516           pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
2517        ENDDO
2518     ENDDO
2519     !cc fin nrlmd le 10/04/2012
2520
2521     if (iflag_thermals>=1) then
2522!jyg<
2523         IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2524!  Appel des thermiques avec les profils exterieurs aux poches
2525          DO k=1,klev
2526           DO i=1,klon
2527            t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2528            q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2529           ENDDO
2530          ENDDO
2531         ELSE
2532!  Appel des thermiques avec les profils moyens
2533          DO k=1,klev
2534           DO i=1,klon
2535            t_therm(i,k) = t_seri(i,k)
2536            q_therm(i,k) = q_seri(i,k)
2537           ENDDO
2538          ENDDO
2539         ENDIF
2540!>jyg
2541        call calltherm(pdtphys &
2542             ,pplay,paprs,pphi,weak_inversion &
2543!!             ,u_seri,v_seri,t_seri,q_seri,zqsat,debut &  !jyg
2544             ,u_seri,v_seri,t_therm,q_therm,zqsat,debut &  !jyg
2545             ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
2546             ,fm_therm,entr_therm,detr_therm &
2547             ,zqasc,clwcon0th,lmax_th,ratqscth &
2548             ,ratqsdiff,zqsatth &
2549             !on rajoute ale et alp, et les caracteristiques de la couche alim
2550             ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
2551             ,ztv,zpspsk,ztla,zthl &
2552             !cc nrlmd le 10/04/2012
2553             ,pbl_tke_input,pctsrf,omega,airephy &
2554             ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
2555             ,n2,s2,ale_bl_stat &
2556             ,therm_tke_max,env_tke_max &
2557             ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
2558             ,alp_bl_conv,alp_bl_stat &
2559             !cc fin nrlmd le 10/04/2012
2560             ,zqla,ztva )
2561!
2562!jyg<
2563         IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2564!  Si les thermiques ne sont presents que hors des poches, la tendance moyenne
2565!  associée doit etre multipliee par la fraction surfacique qu'ils couvrent.
2566          DO k=1,klev
2567           DO i=1,klon
2568!
2569            wake_deltat(i,k) = wake_deltat(i,k) - d_t_ajs(i,k)
2570            wake_deltaq(i,k) = wake_deltaq(i,k) - d_q_ajs(i,k)
2571            t_seri(i,k) = t_therm(i,k) + wake_s(i)*wake_deltat(i,k)
2572            q_seri(i,k) = q_therm(i,k) + wake_s(i)*wake_deltaq(i,k)
2573!
2574            d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
2575            d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
2576            d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
2577            d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
2578!
2579           ENDDO
2580          ENDDO
2581         ELSE
2582          DO k=1,klev
2583           DO i=1,klon
2584            t_seri(i,k) = t_therm(i,k)
2585            q_seri(i,k) = q_therm(i,k)
2586           ENDDO
2587          ENDDO
2588         ENDIF
2589!>jyg
2590
2591        !cc nrlmd le 10/04/2012
2592        !-----------Stochastic triggering-----------
2593        if (iflag_trig_bl.ge.1) then
2594           !
2595           IF (prt_level .GE. 10) THEN
2596              print *,'cin, ale_bl_stat, alp_bl_stat ', &
2597                   cin, ale_bl_stat, alp_bl_stat
2598           ENDIF
2599
2600
2601           !----Initialisations
2602           do i=1,klon
2603              proba_notrig(i)=1.
2604              random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i))
2605              if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
2606                 tau_trig(i)=tau_trig_shallow
2607              else
2608                 tau_trig(i)=tau_trig_deep
2609              endif
2610           enddo
2611           !
2612           IF (prt_level .GE. 10) THEN
2613              print *,'random_notrig, tau_trig ', &
2614                   random_notrig, tau_trig
2615              print *,'s_trig,s2,n2 ', &
2616                   s_trig,s2,n2
2617           ENDIF
2618
2619           !Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
2620           IF (iflag_trig_bl.eq.1) then
2621
2622              !----Tirage al\'eatoire et calcul de ale_bl_trig
2623              do i=1,klon
2624                 if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) )  then
2625                    proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2626                         (n2(i)*dtime/tau_trig(i))
2627                    !        print *, 'proba_notrig(i) ',proba_notrig(i)
2628                    if (random_notrig(i) .ge. proba_notrig(i)) then
2629                       ale_bl_trig(i)=ale_bl_stat(i)
2630                    else
2631                       ale_bl_trig(i)=0.
2632                    endif
2633                 else
2634                    proba_notrig(i)=1.
2635                    random_notrig(i)=0.
2636                    ale_bl_trig(i)=0.
2637                 endif
2638              enddo
2639
2640           ELSE IF (iflag_trig_bl.ge.2) then
2641
2642              do i=1,klon
2643                 if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) )  then
2644                    proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2645                         (n2(i)*dtime/tau_trig(i))
2646                    !        print *, 'proba_notrig(i) ',proba_notrig(i)
2647                    if (random_notrig(i) .ge. proba_notrig(i)) then
2648                       ale_bl_trig(i)=Ale_bl(i)
2649                    else
2650                       ale_bl_trig(i)=0.
2651                    endif
2652                 else
2653                    proba_notrig(i)=1.
2654                    random_notrig(i)=0.
2655                    ale_bl_trig(i)=0.
2656                 endif
2657              enddo
2658
2659           ENDIF
2660
2661           !
2662           IF (prt_level .GE. 10) THEN
2663              print *,'proba_notrig, ale_bl_trig ', &
2664                   proba_notrig, ale_bl_trig
2665           ENDIF
2666
2667        endif !(iflag_trig_bl)
2668
2669        !-----------Statistical closure-----------
2670        if (iflag_clos_bl.eq.1) then
2671
2672           do i=1,klon
2673              !CR: alp probabiliste
2674              if (ale_bl_trig(i).gt.0.) then
2675                 alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
2676              endif
2677           enddo
2678
2679        else if (iflag_clos_bl.eq.2) then
2680
2681           !CR: alp calculee dans thermcell_main
2682           do i=1,klon
2683              alp_bl(i)=alp_bl_stat(i)
2684           enddo
2685
2686        else
2687
2688           alp_bl_stat(:)=0.
2689
2690        endif !(iflag_clos_bl)
2691
2692        IF (prt_level .GE. 10) THEN
2693           print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
2694        ENDIF
2695
2696        !cc fin nrlmd le 10/04/2012
2697
2698        ! ----------------------------------------------------------------------
2699        ! Transport de la TKE par les panaches thermiques.
2700        ! FH : 2010/02/01
2701        !     if (iflag_pbl.eq.10) then
2702        !     call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
2703        !    s           rg,paprs,pbl_tke)
2704        !     endif
2705        ! ----------------------------------------------------------------------
2706        !IM/FH: 2011/02/23
2707        ! Couplage Thermiques/Emanuel seulement si T<0
2708        if (iflag_coupl==2) then
2709         IF (prt_level .GE. 10) THEN
2710           print*,'Couplage Thermiques/Emanuel seulement si T<0'
2711         ENDIF
2712           do i=1,klon
2713              if (t_seri(i,lmax_th(i))>273.) then
2714                 Ale_bl(i)=0.
2715              endif
2716           enddo
2717        endif
2718
2719        do i=1,klon
2720           !           zmax_th(i)=pphi(i,lmax_th(i))/rg
2721           !CR:04/05/12:correction calcul zmax
2722           zmax_th(i)=zmax0(i)
2723        enddo
2724
2725     endif
2726
2727
2728     !  Ajustement sec
2729     !  ==============
2730
2731     ! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
2732     ! a partir du sommet des thermiques.
2733     ! Dans le cas contraire, on demarre au niveau 1.
2734
2735     if (iflag_thermals>=13.or.iflag_thermals<=0) then
2736
2737        if(iflag_thermals.eq.0) then
2738           IF(prt_level>9)WRITE(lunout,*)'ajsec'
2739           limbas(:)=1
2740        else
2741           limbas(:)=lmax_th(:)
2742        endif
2743
2744        ! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2745        ! pour des test de convergence numerique.
2746        ! Le nouveau ajsec est a priori mieux, meme pour le cas
2747        ! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2748        ! non nulles numeriquement pour des mailles non concernees.
2749
2750        if (iflag_thermals==0) then
2751           ! Calling adjustment alone (but not the thermal plume model)
2752           CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
2753                , d_t_ajsb, d_q_ajsb)
2754        else if (iflag_thermals>0) then
2755           ! Calling adjustment above the top of thermal plumes
2756           CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
2757                , d_t_ajsb, d_q_ajsb)
2758        endif
2759
2760        !-----------------------------------------------------------------------
2761        ! ajout des tendances de l'ajustement sec ou des thermiques
2762        CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs,'ajsb')
2763        d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2764        d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2765
2766        !---------------------------------------------------------------------
2767
2768     endif
2769
2770  endif
2771  !
2772  !===================================================================
2773  !IM
2774  IF (ip_ebil_phy.ge.2) THEN
2775     ztit='after dry_adjust'
2776     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2777          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2778          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2779     call diagphy(airephy,ztit,ip_ebil_phy &
2780          , zero_v, zero_v, zero_v, zero_v, zero_v &
2781          , zero_v, zero_v, zero_v, ztsol &
2782          , d_h_vcol, d_qt, d_ec &
2783          , fs_bound, fq_bound )
2784  END IF
2785
2786
2787  !-------------------------------------------------------------------------
2788  ! Computation of ratqs, the width (normalized) of the subrid scale
2789  ! water distribution
2790  CALL  calcratqs(klon,klev,prt_level,lunout,        &
2791       iflag_ratqs,iflag_con,iflag_cldth,pdtphys,  &
2792       ratqsbas,ratqshaut,tau_ratqs,fact_cldcon,   &
2793       ptconv,ptconvth,clwcon0th, rnebcon0th,     &
2794       paprs,pplay,q_seri,zqsat,fm_therm, &
2795       ratqs,ratqsc)
2796
2797
2798  !
2799  ! Appeler le processus de condensation a grande echelle
2800  ! et le processus de precipitation
2801  !-------------------------------------------------------------------------
2802  IF (prt_level .GE.10) THEN
2803     print *,'itap, ->fisrtilp ',itap
2804  ENDIF
2805  !
2806  CALL fisrtilp(dtime,paprs,pplay, &
2807       t_seri, q_seri,ptconv,ratqs, &
2808       d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
2809       rain_lsc, snow_lsc, &
2810       pfrac_impa, pfrac_nucl, pfrac_1nucl, &
2811       frac_impa, frac_nucl, beta_prec_fisrt, &
2812       prfl, psfl, rhcl,  &
2813       zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldth, &
2814       iflag_ice_thermo)
2815  !
2816  WHERE (rain_lsc < 0) rain_lsc = 0.
2817  WHERE (snow_lsc < 0) snow_lsc = 0.
2818
2819  CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs,'lsc')
2820  !---------------------------------------------------------------------------
2821  DO k = 1, klev
2822     DO i = 1, klon
2823        cldfra(i,k) = rneb(i,k)
2824!CR: a quoi ca sert? Faut-il ajouter qs_seri?
2825        IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2826     ENDDO
2827  ENDDO
2828  IF (check) THEN
2829     za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2830     WRITE(lunout,*)"apresilp=", za
2831     zx_t = 0.0
2832     za = 0.0
2833     DO i = 1, klon
2834        za = za + airephy(i)/REAL(klon)
2835        zx_t = zx_t + (rain_lsc(i) &
2836             + snow_lsc(i))*airephy(i)/REAL(klon)
2837     ENDDO
2838     zx_t = zx_t/za*dtime
2839     WRITE(lunout,*)"Precip=", zx_t
2840  ENDIF
2841  !IM
2842  IF (ip_ebil_phy.ge.2) THEN
2843     ztit='after fisrt'
2844     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2845          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2846          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2847     call diagphy(airephy,ztit,ip_ebil_phy &
2848          , zero_v, zero_v, zero_v, zero_v, zero_v &
2849          , zero_v, rain_lsc, snow_lsc, ztsol &
2850          , d_h_vcol, d_qt, d_ec &
2851          , fs_bound, fq_bound )
2852  END IF
2853
2854  if (mydebug) then
2855     call writefield_phy('u_seri',u_seri,llm)
2856     call writefield_phy('v_seri',v_seri,llm)
2857     call writefield_phy('t_seri',t_seri,llm)
2858     call writefield_phy('q_seri',q_seri,llm)
2859  endif
2860
2861  !
2862  !-------------------------------------------------------------------
2863  !  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2864  !-------------------------------------------------------------------
2865
2866  ! 1. NUAGES CONVECTIFS
2867  !
2868  !IM cf FH
2869  !     IF (iflag_cldth.eq.-1) THEN ! seulement pour Tiedtke
2870  IF (iflag_cldth.le.-1) THEN ! seulement pour Tiedtke
2871     snow_tiedtke=0.
2872     !     print*,'avant calcul de la pseudo precip '
2873     !     print*,'iflag_cldth',iflag_cldth
2874     if (iflag_cldth.eq.-1) then
2875        rain_tiedtke=rain_con
2876     else
2877        !       print*,'calcul de la pseudo precip '
2878        rain_tiedtke=0.
2879        !         print*,'calcul de la pseudo precip 0'
2880        do k=1,klev
2881           do i=1,klon
2882              if (d_q_con(i,k).lt.0.) then
2883                 rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
2884                      *(paprs(i,k)-paprs(i,k+1))/rg
2885              endif
2886           enddo
2887        enddo
2888     endif
2889     !
2890     !     call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2891     !
2892
2893     ! Nuages diagnostiques pour Tiedtke
2894     CALL diagcld1(paprs,pplay, &
2895          !IM cf FH  .             rain_con,snow_con,ibas_con,itop_con,
2896          rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
2897          diafra,dialiq)
2898     DO k = 1, klev
2899        DO i = 1, klon
2900           IF (diafra(i,k).GT.cldfra(i,k)) THEN
2901              cldliq(i,k) = dialiq(i,k)
2902              cldfra(i,k) = diafra(i,k)
2903           ENDIF
2904        ENDDO
2905     ENDDO
2906
2907  ELSE IF (iflag_cldth.ge.3) THEN
2908     !  On prend pour les nuages convectifs le max du calcul de la
2909     !  convection et du calcul du pas de temps precedent diminue d'un facteur
2910     !  facttemps
2911     facteur = pdtphys *facttemps
2912     do k=1,klev
2913        do i=1,klon
2914           rnebcon(i,k)=rnebcon(i,k)*facteur
2915           if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k)) &
2916                then
2917              rnebcon(i,k)=rnebcon0(i,k)
2918              clwcon(i,k)=clwcon0(i,k)
2919           endif
2920        enddo
2921     enddo
2922
2923     !
2924     !jq - introduce the aerosol direct and first indirect radiative forcings
2925     !jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2926     IF (flag_aerosol .gt. 0) THEN
2927        IF (iflag_rrtm .EQ. 0) THEN !--old radiation
2928           IF (.NOT. aerosol_couple) THEN
2929              !
2930              CALL readaerosol_optic( &
2931                   debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
2932                   pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
2933                   mass_solu_aero, mass_solu_aero_pi,  &
2934                   tau_aero, piz_aero, cg_aero,  &
2935                   tausum_aero, tau3d_aero)
2936           ENDIF
2937         ELSE                       ! RRTM radiation
2938           IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
2939            abort_message='config_inca=aero et rrtm=1 impossible'
2940            call abort_gcm(modname,abort_message,1)
2941           ELSE
2942!
2943#ifdef CPP_RRTM
2944             CALL readaerosol_optic_rrtm( debut, aerosol_couple, &
2945             new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
2946             pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
2947             tr_seri, mass_solu_aero, mass_solu_aero_pi,  &
2948             tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,  &
2949             tausum_aero, tau3d_aero)
2950
2951             CALL aeropt_lw_rrtm
2952#else
2953
2954              abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
2955              call abort_gcm(modname,abort_message,1)
2956#endif
2957              !
2958           ENDIF
2959        ENDIF
2960     ELSE
2961        tausum_aero(:,:,:) = 0.
2962        IF (iflag_rrtm .EQ. 0) THEN !--old radiation
2963           tau_aero(:,:,:,:) = 0.
2964           piz_aero(:,:,:,:) = 0.
2965           cg_aero(:,:,:,:)  = 0.
2966        ELSE
2967           tau_aero_sw_rrtm(:,:,:,:)=0.0
2968           piz_aero_sw_rrtm(:,:,:,:)=0.0
2969           cg_aero_sw_rrtm(:,:,:,:)=0.0
2970        ENDIF
2971     ENDIF
2972     !
2973     !--STRAT AEROSOL
2974     !--updates tausum_aero,tau_aero,piz_aero,cg_aero
2975     IF (flag_aerosol_strat) THEN
2976        IF (prt_level .GE.10) THEN
2977         PRINT *,'appel a readaerosolstrat', mth_cur
2978        ENDIF
2979        IF (iflag_rrtm.EQ.0) THEN
2980           CALL readaerosolstrato(debut)
2981        ELSE
2982#ifdef CPP_RRTM
2983           CALL readaerosolstrato_rrtm(debut)
2984#else
2985
2986           abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
2987           call abort_gcm(modname,abort_message,1)
2988#endif
2989        ENDIF
2990     ENDIF
2991     !--fin STRAT AEROSOL
2992
2993
2994     !   On prend la somme des fractions nuageuses et des contenus en eau
2995
2996     if (iflag_cldth>=5) then
2997
2998        do k=1,klev
2999           ptconvth(:,k)=fm_therm(:,k+1)>0.
3000        enddo
3001
3002        if (iflag_coupl==4) then
3003
3004           ! Dans le cas iflag_coupl==4, on prend la somme des convertures
3005           ! convectives et lsc dans la partie des thermiques
3006           ! Le controle par iflag_coupl est peut etre provisoire.
3007           do k=1,klev
3008              do i=1,klon
3009                 if (ptconv(i,k).and.ptconvth(i,k)) then
3010                    cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3011                    cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3012                 else if (ptconv(i,k)) then
3013                    cldfra(i,k)=rnebcon(i,k)
3014                    cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3015                 endif
3016              enddo
3017           enddo
3018
3019        else if (iflag_coupl==5) then
3020           do k=1,klev
3021              do i=1,klon
3022                 cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3023                 cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3024              enddo
3025           enddo
3026
3027        else
3028
3029           ! Si on est sur un point touche par la convection profonde et pas
3030           ! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
3031           ! de la convection profonde.
3032
3033           !IM/FH: 2011/02/23
3034           ! definition des points sur lesquels ls thermiques sont actifs
3035
3036           do k=1,klev
3037              do i=1,klon
3038                 if (ptconv(i,k).and. .not. ptconvth(i,k)) then
3039                    cldfra(i,k)=rnebcon(i,k)
3040                    cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3041                 endif
3042              enddo
3043           enddo
3044
3045        endif
3046
3047     else
3048
3049        ! Ancienne version
3050        cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3051        cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3052     endif
3053
3054  ENDIF
3055
3056  !     plulsc(:)=0.
3057  !     do k=1,klev,-1
3058  !        do i=1,klon
3059  !              zzz=prfl(:,k)+psfl(:,k)
3060  !           if (.not.ptconvth.zzz.gt.0.)
3061  !        enddo prfl, psfl,
3062  !     enddo
3063  !
3064  ! 2. NUAGES STARTIFORMES
3065  !
3066  IF (ok_stratus) THEN
3067     CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3068     DO k = 1, klev
3069        DO i = 1, klon
3070           IF (diafra(i,k).GT.cldfra(i,k)) THEN
3071              cldliq(i,k) = dialiq(i,k)
3072              cldfra(i,k) = diafra(i,k)
3073           ENDIF
3074        ENDDO
3075     ENDDO
3076  ENDIF
3077  !
3078  ! Precipitation totale
3079  !
3080  DO i = 1, klon
3081     rain_fall(i) = rain_con(i) + rain_lsc(i)
3082     snow_fall(i) = snow_con(i) + snow_lsc(i)
3083  ENDDO
3084  !IM
3085  IF (ip_ebil_phy.ge.2) THEN
3086     ztit="after diagcld"
3087     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3088          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3089          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3090     call diagphy(airephy,ztit,ip_ebil_phy &
3091          , zero_v, zero_v, zero_v, zero_v, zero_v &
3092          , zero_v, zero_v, zero_v, ztsol &
3093          , d_h_vcol, d_qt, d_ec &
3094          , fs_bound, fq_bound )
3095  END IF
3096  !
3097  ! Calculer l'humidite relative pour diagnostique
3098  !
3099  DO k = 1, klev
3100     DO i = 1, klon
3101        zx_t = t_seri(i,k)
3102        IF (thermcep) THEN
3103           if (iflag_ice_thermo.eq.0) then
3104           zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3105           else
3106           zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t))
3107           endif
3108           zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3109           zx_qs  = MIN(0.5,zx_qs)
3110           zcor   = 1./(1.-retv*zx_qs)
3111           zx_qs  = zx_qs*zcor
3112        ELSE
3113           IF (zx_t.LT.t_coup) THEN
3114              zx_qs = qsats(zx_t)/pplay(i,k)
3115           ELSE
3116              zx_qs = qsatl(zx_t)/pplay(i,k)
3117           ENDIF
3118        ENDIF
3119        zx_rh(i,k) = q_seri(i,k)/zx_qs
3120        zqsat(i,k)=zx_qs
3121     ENDDO
3122  ENDDO
3123
3124  !IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3125  !   equivalente a 2m (tpote) pour diagnostique
3126  !
3127  DO i = 1, klon
3128     tpot(i)=zt2m(i)*(100000./paprs(i,1))**RKAPPA
3129     IF (thermcep) THEN
3130        IF(zt2m(i).LT.RTT) then
3131           Lheat=RLSTT
3132        ELSE
3133           Lheat=RLVTT
3134        ENDIF
3135     ELSE
3136        IF (zt2m(i).LT.RTT) THEN
3137           Lheat=RLSTT
3138        ELSE
3139           Lheat=RLVTT
3140        ENDIF
3141     ENDIF
3142     tpote(i) = tpot(i)*      &
3143          EXP((Lheat *qsat2m(i))/(RCPD*zt2m(i)))
3144  ENDDO
3145
3146  IF (type_trac == 'inca') THEN
3147#ifdef INCA
3148     CALL VTe(VTphysiq)
3149     CALL VTb(VTinca)
3150     calday = REAL(days_elapsed + 1) + jH_cur
3151
3152     call chemtime(itap+itau_phy-1, date0, dtime, itap)
3153     IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
3154        CALL AEROSOL_METEO_CALC( &
3155             calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3156             prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
3157     END IF
3158
3159     zxsnow_dummy(:) = 0.0
3160
3161     CALL chemhook_begin (calday, &
3162          days_elapsed+1, &
3163          jH_cur, &
3164          pctsrf(1,1), &
3165          rlat, &
3166          rlon, &
3167          airephy, &
3168          paprs, &
3169          pplay, &
3170          coefh(1:klon,1:klev,is_ave), &
3171          pphi, &
3172          t_seri, &
3173          u, &
3174          v, &
3175          wo(:, :, 1), &
3176          q_seri, &
3177          zxtsol, &
3178          zxsnow_dummy, &
3179          solsw, &
3180          albsol1, &
3181          rain_fall, &
3182          snow_fall, &
3183          itop_con, &
3184          ibas_con, &
3185          cldfra, &
3186          iim, &
3187          jjm, &
3188          tr_seri, &
3189          ftsol, &
3190          paprs, &
3191          cdragh, &
3192          cdragm, &
3193          pctsrf, &
3194          pdtphys, &
3195          itap)
3196
3197     CALL VTe(VTinca)
3198     CALL VTb(VTphysiq)
3199#endif
3200  END IF !type_trac = inca
3201  !     
3202  ! Calculer les parametres optiques des nuages et quelques
3203  ! parametres pour diagnostiques:
3204  !
3205
3206  IF (aerosol_couple) THEN
3207     mass_solu_aero(:,:)    = ccm(:,:,1)
3208     mass_solu_aero_pi(:,:) = ccm(:,:,2)
3209  END IF
3210
3211  if (ok_newmicro) then
3212     IF (iflag_rrtm.NE.0) THEN
3213#ifdef CPP_RRTM
3214        IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3215           abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 pour ok_cdnc'
3216           call abort_gcm(modname,abort_message,1)
3217        endif
3218#else
3219
3220        abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
3221        call abort_gcm(modname,abort_message,1)
3222#endif
3223     ENDIF
3224     CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
3225          paprs, pplay, t_seri, cldliq, cldfra, &
3226          cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3227          flwp, fiwp, flwc, fiwc, &
3228          mass_solu_aero, mass_solu_aero_pi, &
3229          cldtaupi, re, fl, ref_liq, ref_ice, &
3230          ref_liq_pi, ref_ice_pi)
3231  else
3232     CALL nuage (paprs, pplay, &
3233          t_seri, cldliq, cldfra, cldtau, cldemi, &
3234          cldh, cldl, cldm, cldt, cldq, &
3235          ok_aie, &
3236          mass_solu_aero, mass_solu_aero_pi, &
3237          bl95_b0, bl95_b1, &
3238          cldtaupi, re, fl)
3239  endif
3240  !
3241  !IM betaCRF
3242  !
3243  cldtaurad   = cldtau
3244  cldtaupirad = cldtaupi
3245  cldemirad   = cldemi
3246
3247  !
3248  if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3249       lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3250     !
3251     ! global
3252     !
3253     DO k=1, klev
3254        DO i=1, klon
3255           if (pplay(i,k).GE.pfree) THEN
3256              beta(i,k) = beta_pbl
3257           else
3258              beta(i,k) = beta_free
3259           endif
3260           if (mskocean_beta) THEN
3261              beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3262           endif
3263           cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
3264           cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3265           cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
3266           cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
3267        ENDDO
3268     ENDDO
3269     !
3270  else
3271     !
3272     ! regional
3273     !
3274     DO k=1, klev
3275        DO i=1,klon
3276           !
3277           if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND. &
3278                rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
3279              if (pplay(i,k).GE.pfree) THEN
3280                 beta(i,k) = beta_pbl
3281              else
3282                 beta(i,k) = beta_free
3283              endif
3284              if (mskocean_beta) THEN
3285                 beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3286              endif
3287              cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
3288              cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3289              cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
3290              cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
3291           endif
3292           !
3293        ENDDO
3294     ENDDO
3295     !
3296  endif
3297  !
3298  ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3299  !
3300  IF (MOD(itaprad,radpas).EQ.0) THEN
3301
3302     DO i = 1, klon
3303        albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce) &
3304             + falb1(i,is_lic) * pctsrf(i,is_lic) &
3305             + falb1(i,is_ter) * pctsrf(i,is_ter) &
3306             + falb1(i,is_sic) * pctsrf(i,is_sic)
3307        albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce) &
3308             + falb2(i,is_lic) * pctsrf(i,is_lic) &
3309             + falb2(i,is_ter) * pctsrf(i,is_ter) &
3310             + falb2(i,is_sic) * pctsrf(i,is_sic)
3311     ENDDO
3312
3313     if (mydebug) then
3314        call writefield_phy('u_seri',u_seri,llm)
3315        call writefield_phy('v_seri',v_seri,llm)
3316        call writefield_phy('t_seri',t_seri,llm)
3317        call writefield_phy('q_seri',q_seri,llm)
3318     endif
3319
3320     IF (aerosol_couple) THEN
3321#ifdef INCA
3322        CALL radlwsw_inca  &
3323             (kdlon,kflev,dist, rmu0, fract, solaire, &
3324             paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3325             wo(:, :, 1), &
3326             cldfrarad, cldemirad, cldtaurad, &
3327             heat,heat0,cool,cool0,radsol,albpla, &
3328             topsw,toplw,solsw,sollw, &
3329             sollwdown, &
3330             topsw0,toplw0,solsw0,sollw0, &
3331             lwdn0, lwdn, lwup0, lwup,  &
3332             swdn0, swdn, swup0, swup, &
3333             ok_ade, ok_aie, &
3334             tau_aero, piz_aero, cg_aero, &
3335             topswad_aero, solswad_aero, &
3336             topswad0_aero, solswad0_aero, &
3337             topsw_aero, topsw0_aero, &
3338             solsw_aero, solsw0_aero, &
3339             cldtaupirad, &
3340             topswai_aero, solswai_aero)
3341
3342#endif
3343     ELSE
3344        !
3345        !IM calcul radiatif pour le cas actuel
3346        !
3347        RCO2 = RCO2_act
3348        RCH4 = RCH4_act
3349        RN2O = RN2O_act
3350        RCFC11 = RCFC11_act
3351        RCFC12 = RCFC12_act
3352        !
3353        IF (prt_level .GE.10) THEN
3354           print *,' ->radlwsw, number 1 '
3355        ENDIF
3356        !
3357        CALL radlwsw &
3358             (dist, rmu0, fract,  &
3359             paprs, pplay,zxtsol,albsol1, albsol2,  &
3360             t_seri,q_seri,wo, &
3361             cldfrarad, cldemirad, cldtaurad, &
3362             ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3363             flag_aerosol_strat, &
3364             tau_aero, piz_aero, cg_aero, &
3365             tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,&     ! Rajoute par OB pour RRTM
3366             tau_aero_lw_rrtm, &
3367             cldtaupirad,new_aod, &
3368             zqsat, flwc, fiwc, &
3369             ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3370             heat,heat0,cool,cool0,radsol,albpla, &
3371             topsw,toplw,solsw,sollw, &
3372             sollwdown, &
3373             topsw0,toplw0,solsw0,sollw0, &
3374             lwdn0, lwdn, lwup0, lwup,  &
3375             swdn0, swdn, swup0, swup, &
3376             topswad_aero, solswad_aero, &
3377             topswai_aero, solswai_aero, &
3378             topswad0_aero, solswad0_aero, &
3379             topsw_aero, topsw0_aero, &
3380             solsw_aero, solsw0_aero, &
3381             topswcf_aero, solswcf_aero, &
3382             !-C. Kleinschmitt for LW diagnostics
3383             toplwad_aero, sollwad_aero,&
3384             toplwai_aero, sollwai_aero, &
3385             toplwad0_aero, sollwad0_aero,&
3386             !-end
3387             ZLWFT0_i, ZFLDN0, ZFLUP0, &
3388             ZSWFT0_i, ZFSDN0, ZFSUP0)
3389
3390        !
3391        !IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3392        !IM des taux doit etre different du taux actuel
3393        !IM Par defaut on a les taux perturbes egaux aux taux actuels
3394        !
3395        if (ok_4xCO2atm) then
3396           if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
3397                RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
3398                RCFC12_per.NE.RCFC12_act) THEN
3399              !
3400              RCO2 = RCO2_per
3401              RCH4 = RCH4_per
3402              RN2O = RN2O_per
3403              RCFC11 = RCFC11_per
3404              RCFC12 = RCFC12_per
3405              !
3406              IF (prt_level .GE.10) THEN
3407                 print *,' ->radlwsw, number 2 '
3408              ENDIF
3409              !
3410              CALL radlwsw &
3411                   (dist, rmu0, fract,  &
3412                   paprs, pplay,zxtsol,albsol1, albsol2,  &
3413                   t_seri,q_seri,wo, &
3414                   cldfra, cldemi, cldtau, &
3415                   ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3416                   flag_aerosol_strat, &
3417                   tau_aero, piz_aero, cg_aero, &
3418                   tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,&     ! Rajoute par OB pour RRTM
3419                   tau_aero_lw_rrtm, &
3420                   cldtaupi,new_aod, &
3421                   zqsat, flwc, fiwc, &
3422                   ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3423                   heatp,heat0p,coolp,cool0p,radsolp,albplap, &
3424                   topswp,toplwp,solswp,sollwp, &
3425                   sollwdownp, &
3426                   topsw0p,toplw0p,solsw0p,sollw0p, &
3427                   lwdn0p, lwdnp, lwup0p, lwupp,  &
3428                   swdn0p, swdnp, swup0p, swupp, &
3429                   topswad_aerop, solswad_aerop, &
3430                   topswai_aerop, solswai_aerop, &
3431                   topswad0_aerop, solswad0_aerop, &
3432                   topsw_aerop, topsw0_aerop, &
3433                   solsw_aerop, solsw0_aerop, &
3434                   topswcf_aerop, solswcf_aerop, &
3435                   !-C. Kleinschmitt for LW diagnostics
3436                   toplwad_aerop, sollwad_aerop,&
3437                   toplwai_aerop, sollwai_aerop, &
3438                   toplwad0_aerop, sollwad0_aerop,&
3439                   !-end
3440                   ZLWFT0_i, ZFLDN0, ZFLUP0, &
3441                   ZSWFT0_i, ZFSDN0, ZFSUP0)
3442           endif
3443        endif
3444        !
3445     ENDIF ! aerosol_couple
3446     itaprad = 0
3447  ENDIF ! MOD(itaprad,radpas)
3448  itaprad = itaprad + 1
3449
3450  IF (iflag_radia.eq.0) THEN
3451     IF (prt_level.ge.9) THEN
3452        PRINT *,'--------------------------------------------------'
3453        PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3454        PRINT *,'>>>>           heat et cool mis a zero '
3455        PRINT *,'--------------------------------------------------'
3456     END IF
3457     heat=0.
3458     cool=0.
3459     sollw=0.   ! MPL 01032011
3460     solsw=0.
3461     radsol=0.
3462     swup=0.    ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
3463     swup0=0.
3464     swdn=0.
3465     swdn0=0.
3466     lwup=0.
3467     lwup0=0.
3468     lwdn=0.
3469     lwdn0=0.
3470  END IF
3471
3472  !
3473  ! Ajouter la tendance des rayonnements (tous les pas)
3474  !
3475  d_t_swr(:,:)=heat(:,:)*dtime/RDAY
3476  d_t_lwr(:,:)=-cool(:,:)*dtime/RDAY
3477  d_t_sw0(:,:)=heat0(:,:)*dtime/RDAY
3478  d_t_lw0(:,:)=-cool0(:,:)*dtime/RDAY
3479  CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW')
3480  CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW')
3481
3482  !
3483  if (mydebug) then
3484     call writefield_phy('u_seri',u_seri,llm)
3485     call writefield_phy('v_seri',v_seri,llm)
3486     call writefield_phy('t_seri',t_seri,llm)
3487     call writefield_phy('q_seri',q_seri,llm)
3488  endif
3489
3490  !IM
3491  IF (ip_ebil_phy.ge.2) THEN
3492     ztit='after rad'
3493     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3494          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3495          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3496     call diagphy(airephy,ztit,ip_ebil_phy &
3497          , topsw, toplw, solsw, sollw, zero_v &
3498          , zero_v, zero_v, zero_v, ztsol &
3499          , d_h_vcol, d_qt, d_ec &
3500          , fs_bound, fq_bound )
3501  END IF
3502  !
3503  !
3504  ! Calculer l'hydrologie de la surface
3505  !
3506  !      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3507  !     .            agesno, ftsol,fqsurf,fsnow, ruis)
3508  !
3509
3510  !
3511  ! Calculer le bilan du sol et la derive de temperature (couplage)
3512  !
3513  DO i = 1, klon
3514     !         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3515     ! a la demande de JLD
3516     bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3517  ENDDO
3518  !
3519  !moddeblott(jan95)
3520  ! Appeler le programme de parametrisation de l'orographie
3521  ! a l'echelle sous-maille:
3522  !
3523  IF (prt_level .GE.10) THEN
3524     print *,' call orography ? ', ok_orodr
3525  ENDIF
3526  !
3527  IF (ok_orodr) THEN
3528     !
3529     !  selection des points pour lesquels le shema est actif:
3530     igwd=0
3531     DO i=1,klon
3532        itest(i)=0
3533        !        IF ((zstd(i).gt.10.0)) THEN
3534        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3535           itest(i)=1
3536           igwd=igwd+1
3537           idx(igwd)=i
3538        ENDIF
3539     ENDDO
3540     !        igwdim=MAX(1,igwd)
3541     !
3542     IF (ok_strato) THEN
3543
3544        CALL drag_noro_strato(klon,klev,dtime,paprs,pplay, &
3545             zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3546             igwd,idx,itest, &
3547             t_seri, u_seri, v_seri, &
3548             zulow, zvlow, zustrdr, zvstrdr, &
3549             d_t_oro, d_u_oro, d_v_oro)
3550
3551     ELSE
3552        CALL drag_noro(klon,klev,dtime,paprs,pplay, &
3553             zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3554             igwd,idx,itest, &
3555             t_seri, u_seri, v_seri, &
3556             zulow, zvlow, zustrdr, zvstrdr, &
3557             d_t_oro, d_u_oro, d_v_oro)
3558     ENDIF
3559     !
3560     !  ajout des tendances
3561     !-----------------------------------------------------------------------------------------
3562     ! ajout des tendances de la trainee de l'orographie
3563     CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro')
3564     !-----------------------------------------------------------------------------------------
3565     !
3566  ENDIF ! fin de test sur ok_orodr
3567  !
3568  if (mydebug) then
3569     call writefield_phy('u_seri',u_seri,llm)
3570     call writefield_phy('v_seri',v_seri,llm)
3571     call writefield_phy('t_seri',t_seri,llm)
3572     call writefield_phy('q_seri',q_seri,llm)
3573  endif
3574
3575  IF (ok_orolf) THEN
3576     !
3577     !  selection des points pour lesquels le shema est actif:
3578     igwd=0
3579     DO i=1,klon
3580        itest(i)=0
3581        IF ((zpic(i)-zmea(i)).GT.100.) THEN
3582           itest(i)=1
3583           igwd=igwd+1
3584           idx(igwd)=i
3585        ENDIF
3586     ENDDO
3587     !        igwdim=MAX(1,igwd)
3588     !
3589     IF (ok_strato) THEN
3590
3591        CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
3592             rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
3593             igwd,idx,itest, &
3594             t_seri, u_seri, v_seri, &
3595             zulow, zvlow, zustrli, zvstrli, &
3596             d_t_lif, d_u_lif, d_v_lif               )
3597
3598     ELSE
3599        CALL lift_noro(klon,klev,dtime,paprs,pplay, &
3600             rlat,zmea,zstd,zpic, &
3601             itest, &
3602             t_seri, u_seri, v_seri, &
3603             zulow, zvlow, zustrli, zvstrli, &
3604             d_t_lif, d_u_lif, d_v_lif)
3605     ENDIF
3606     !   
3607     !-----------------------------------------------------------------------------------------
3608     ! ajout des tendances de la portance de l'orographie
3609     CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,dqi0,paprs,'lif')
3610     !-----------------------------------------------------------------------------------------
3611     !
3612  ENDIF ! fin de test sur ok_orolf
3613  !  HINES GWD PARAMETRIZATION
3614
3615  IF (ok_hines) then
3616
3617     CALL hines_gwd(klon,klev,dtime,paprs,pplay, &
3618          rlat,t_seri,u_seri,v_seri, &
3619          zustrhi,zvstrhi, &
3620          d_t_hin, d_u_hin, d_v_hin)
3621     !
3622     !  ajout des tendances
3623     CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,dqi0,paprs,'hin')
3624
3625  ENDIF
3626
3627  if (ok_gwd_rando) then
3628     call FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
3629          rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
3630          du_gwd_rando, dv_gwd_rando)
3631     CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0,dqi0,paprs, &
3632          'flott_gwd_rando')
3633  end if
3634
3635  ! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3636
3637  if (mydebug) then
3638     call writefield_phy('u_seri',u_seri,llm)
3639     call writefield_phy('v_seri',v_seri,llm)
3640     call writefield_phy('t_seri',t_seri,llm)
3641     call writefield_phy('q_seri',q_seri,llm)
3642  endif
3643
3644  DO i = 1, klon
3645     zustrph(i)=0.
3646     zvstrph(i)=0.
3647  ENDDO
3648  DO k = 1, klev
3649     DO i = 1, klon
3650        zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
3651             (paprs(i,k)-paprs(i,k+1))/rg
3652        zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
3653             (paprs(i,k)-paprs(i,k+1))/rg
3654     ENDDO
3655  ENDDO
3656  !
3657  !IM calcul composantes axiales du moment angulaire et couple des montagnes
3658  !
3659  IF (is_sequential .and. ok_orodr) THEN
3660     CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
3661          ra,rg,romega, &
3662          rlat,rlon,pphis, &
3663          zustrdr,zustrli,zustrph, &
3664          zvstrdr,zvstrli,zvstrph, &
3665          paprs,u,v, &
3666          aam, torsfc)
3667  ENDIF
3668  !IM cf. FLott END
3669  !IM
3670  IF (ip_ebil_phy.ge.2) THEN
3671     ztit='after orography'
3672     CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3673          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3674          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3675     call diagphy(airephy,ztit,ip_ebil_phy &
3676          , zero_v, zero_v, zero_v, zero_v, zero_v &
3677          , zero_v, zero_v, zero_v, ztsol &
3678          , d_h_vcol, d_qt, d_ec &
3679          , fs_bound, fq_bound )
3680  END IF
3681
3682  !DC Calcul de la tendance due au methane
3683  IF(ok_qch4) THEN
3684     CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
3685  ! ajout de la tendance d'humidite due au methane
3686     CALL add_phys_tend(du0,dv0,dt0,d_q_ch4*dtime,dql0,'q_ch4')
3687  END IF
3688  !
3689  !
3690  !====================================================================
3691  ! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3692  !====================================================================
3693  ! Abderrahmane 24.08.09
3694
3695  IF (ok_cosp) THEN
3696     ! adeclarer
3697#ifdef CPP_COSP
3698     IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3699
3700      IF (prt_level .GE.10) THEN
3701        print*,'freq_cosp',freq_cosp
3702      ENDIF
3703        mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3704        !       print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3705        !     s        ref_liq,ref_ice
3706        call phys_cosp(itap,dtime,freq_cosp, &
3707             ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
3708             ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, &
3709             klon,klev,rlon,rlat,presnivs,overlap, &
3710             fract,ref_liq,ref_ice, &
3711             pctsrf(:,is_ter)+pctsrf(:,is_lic), &
3712             zu10m,zv10m,pphis, &
3713             zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
3714             qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
3715             prfl(:,1:klev),psfl(:,1:klev), &
3716             pmflxr(:,1:klev),pmflxs(:,1:klev), &
3717             mr_ozone,cldtau, cldemi)
3718
3719        !     L          calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3720        !     L          cfaddbze,clcalipso2,dbze,cltlidarradar,
3721        !     M          clMISR,
3722        !     R          clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3723        !     I          tauisccp,albisccp,meantbisccp,meantbclrisccp)
3724
3725     ENDIF
3726
3727#endif
3728  ENDIF  !ok_cosp
3729!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3730  !AA
3731  !AA Installation de l'interface online-offline pour traceurs
3732  !AA
3733  !====================================================================
3734  !   Calcul  des tendances traceurs
3735  !====================================================================
3736  !
3737
3738  IF (type_trac=='repr') THEN
3739     sh_in(:,:) = q_seri(:,:)
3740  ELSE
3741     sh_in(:,:) = qx(:,:,ivap)
3742  END IF
3743
3744  call phytrac ( &
3745       itap,     days_elapsed+1,    jH_cur,   debut, &
3746       lafin,    dtime,     u, v,     t, &
3747       paprs,    pplay,     pmfu,     pmfd, &
3748       pen_u,    pde_u,     pen_d,    pde_d, &
3749       cdragh,   coefh(1:klon,1:klev,is_ave),   fm_therm, entr_therm, &
3750       u1,       v1,        ftsol,    pctsrf, &
3751       zustar,   zu10m,     zv10m, &
3752       wstar(:,is_ave),    ale_bl,         ale_wake, &
3753       rlat,     rlon, &
3754       frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
3755       presnivs, pphis,     pphi,     albsol1, &
3756       sh_in,    rhcl,      cldfra,   rneb, &
3757       diafra,   cldliq,    itop_con, ibas_con, &
3758       pmflxr,   pmflxs,    prfl,     psfl, &
3759       da,       phi,       mp,       upwd, &
3760       phi2,     d1a,       dam,      sij, wght_cvfd, &        !<<RomP+RL
3761       wdtrainA, wdtrainM,  sigd,     clw,elij, &   !<<RomP
3762       ev,       ep,        epmlmMm,  eplaMm, &     !<<RomP
3763       dnwd,     aerosol_couple,      flxmass_w, &
3764       tau_aero, piz_aero,  cg_aero,  ccm, &
3765       rfname, &
3766       d_tr_dyn, &                                 !<<RomP
3767       tr_seri)
3768
3769  IF (offline) THEN
3770
3771     IF (prt_level.ge.9) &
3772          print*,'Attention on met a 0 les thermiques pour phystoke'
3773     call phystokenc ( &
3774          nlon,klev,pdtphys,rlon,rlat, &
3775          t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
3776          fm_therm,entr_therm, &
3777          cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
3778          frac_impa, frac_nucl, &
3779          pphis,airephy,dtime,itap, &
3780          qx(:,:,ivap),da,phi,mp,upwd,dnwd)
3781
3782
3783  ENDIF
3784
3785  !
3786  ! Calculer le transport de l'eau et de l'energie (diagnostique)
3787  !
3788  CALL transp (paprs,zxtsol, &
3789       t_seri, q_seri, u_seri, v_seri, zphi, &
3790       ve, vq, ue, uq)
3791  !
3792  !IM global posePB BEG
3793  IF(1.EQ.0) THEN
3794     !
3795     CALL transp_lay (paprs,zxtsol, &
3796          t_seri, q_seri, u_seri, v_seri, zphi, &
3797          ve_lay, vq_lay, ue_lay, uq_lay)
3798     !
3799  ENDIF !(1.EQ.0) THEN
3800  !IM global posePB END
3801  ! Accumuler les variables a stocker dans les fichiers histoire:
3802  !
3803
3804  !================================================================
3805  ! Conversion of kinetic and potential energy into heat, for
3806  ! parameterisation of subgrid-scale motions
3807  !================================================================
3808
3809  d_t_ec(:,:)=0.
3810  forall (k=1: llm) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
3811  CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
3812       u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
3813       zmasse,exner,d_t_ec)
3814  t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
3815
3816  !IM
3817  IF (ip_ebil_phy.ge.1) THEN
3818     ztit='after physic'
3819     CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
3820          , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3821          , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3822     !     Comme les tendances de la physique sont ajoute dans la dynamique,
3823     !     on devrait avoir que la variation d'entalpie par la dynamique
3824     !     est egale a la variation de la physique au pas de temps precedent.
3825     !     Donc la somme de ces 2 variations devrait etre nulle.
3826
3827     call diagphy(airephy,ztit,ip_ebil_phy &
3828          , topsw, toplw, solsw, sollw, sens &
3829          , evap, rain_fall, snow_fall, ztsol &
3830          , d_h_vcol, d_qt, d_ec &
3831          , fs_bound, fq_bound )
3832     !
3833     d_h_vcol_phy=d_h_vcol
3834     !
3835  END IF
3836  !
3837  !=======================================================================
3838  !   SORTIES
3839  !=======================================================================
3840  !
3841  !IM initialisation + calculs divers diag AMIP2
3842  !
3843  include "calcul_divers.h"
3844  !
3845  !IM Interpolation sur les niveaux de pression du NMC
3846  !   -------------------------------------------------
3847  !
3848  include "calcul_STDlev.h"
3849  !
3850  ! slp sea level pressure
3851  slp(:) = paprs(:,1)*exp(pphis(:)/(RD*t_seri(:,1)))
3852  !
3853  !cc prw = eau precipitable
3854  DO i = 1, klon
3855     prw(i) = 0.
3856     DO k = 1, klev
3857        prw(i) = prw(i) + &
3858             q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3859     ENDDO
3860  ENDDO
3861  !
3862  IF (type_trac == 'inca') THEN
3863#ifdef INCA
3864     CALL VTe(VTphysiq)
3865     CALL VTb(VTinca)
3866
3867     CALL chemhook_end ( &
3868          dtime, &
3869          pplay, &
3870          t_seri, &
3871          tr_seri, &
3872          nbtr, &
3873          paprs, &
3874          q_seri, &
3875          airephy, &
3876          pphi, &
3877          pphis, &
3878          zx_rh)
3879
3880     CALL VTe(VTinca)
3881     CALL VTb(VTphysiq)
3882#endif
3883  END IF
3884
3885
3886  !
3887  ! Convertir les incrementations en tendances
3888  !
3889  IF (prt_level .GE.10) THEN
3890     print *,'Convertir les incrementations en tendances '
3891  ENDIF
3892  !
3893  if (mydebug) then
3894     call writefield_phy('u_seri',u_seri,llm)
3895     call writefield_phy('v_seri',v_seri,llm)
3896     call writefield_phy('t_seri',t_seri,llm)
3897     call writefield_phy('q_seri',q_seri,llm)
3898  endif
3899
3900  DO k = 1, klev
3901     DO i = 1, klon
3902        d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3903        d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3904        d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3905        d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3906        d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3907!CR: on ajoute le contenu en glace
3908        if (nqo.eq.3) then
3909        d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
3910        endif
3911     ENDDO
3912  ENDDO
3913  !
3914!CR: nb de traceurs eau: nqo
3915!  IF (nqtot.GE.3) THEN
3916   IF (nqtot.GE.(nqo+1)) THEN
3917!     DO iq = 3, nqtot
3918     DO iq = nqo+1, nqtot
3919        DO  k = 1, klev
3920           DO  i = 1, klon
3921!              d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3922               d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
3923           ENDDO
3924        ENDDO
3925     ENDDO
3926  ENDIF
3927  !
3928  !IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3929  !IM global posePB      include "write_bilKP_ins.h"
3930  !IM global posePB      include "write_bilKP_ave.h"
3931  !
3932
3933  ! Sauvegarder les valeurs de t et q a la fin de la physique:
3934  !
3935  DO k = 1, klev
3936     DO i = 1, klon
3937        u_ancien(i,k) = u_seri(i,k)
3938        v_ancien(i,k) = v_seri(i,k)
3939        t_ancien(i,k) = t_seri(i,k)
3940        q_ancien(i,k) = q_seri(i,k)
3941     ENDDO
3942  ENDDO
3943
3944!!! RomP >>>
3945!CR: nb de traceurs eau: nqo
3946!  IF (nqtot.GE.3) THEN
3947   IF (nqtot.GE.(nqo+1)) THEN
3948!     DO iq = 3, nqtot
3949     DO iq = nqo+1, nqtot
3950        DO k = 1, klev
3951           DO i = 1, klon
3952!              tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
3953              tr_ancien(i,k,iq-nqo) = tr_seri(i,k,iq-nqo)
3954           ENDDO
3955        ENDDO
3956     ENDDO
3957  ENDIF
3958!!! RomP <<<
3959  !==========================================================================
3960  ! Sorties des tendances pour un point particulier
3961  ! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3962  ! pour le debug
3963  ! La valeur de igout est attribuee plus haut dans le programme
3964  !==========================================================================
3965
3966  if (prt_level.ge.1) then
3967     write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3968     write(lunout,*) &
3969          'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
3970     write(lunout,*) &
3971          nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
3972          pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
3973          pctsrf(igout,is_sic)
3974     write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3975     do k=1,klev
3976        write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
3977             d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
3978             d_t_eva(igout,k)
3979     enddo
3980     write(lunout,*) 'cool,heat'
3981     do k=1,klev
3982        write(lunout,*) cool(igout,k),heat(igout,k)
3983     enddo
3984
3985     write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3986     do k=1,klev
3987        write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
3988             d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3989     enddo
3990
3991     write(lunout,*) 'd_ps ',d_ps(igout)
3992     write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3993     do k=1,klev
3994        write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
3995             d_qx(igout,k,1),d_qx(igout,k,2)
3996     enddo
3997  endif
3998
3999  !==========================================================================
4000
4001  !============================================================
4002  !   Calcul de la temperature potentielle
4003  !============================================================
4004  DO k = 1, klev
4005     DO i = 1, klon
4006        !JYG/IM theta en debut du pas de temps
4007        !JYG/IM       theta(i,k)=t(i,k)*(100000./pplay(i,k))**(RD/RCPD)
4008        !JYG/IM theta en fin de pas de temps de physique
4009        theta(i,k)=t_seri(i,k)*(100000./pplay(i,k))**(RD/RCPD)
4010        ! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers     MPL 20130625
4011        ! fth_fonctions.F90 et parkind1.F90
4012        ! sinon thetal=theta
4013        !       thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4014        !    :         ql_seri(i,k))
4015        thetal(i,k)=theta(i,k)
4016     ENDDO
4017  ENDDO
4018  !
4019
4020  ! 22.03.04 BEG
4021  !=============================================================
4022  !   Ecriture des sorties
4023  !=============================================================
4024#ifdef CPP_IOIPSL
4025
4026  ! Recupere des varibles calcule dans differents modules
4027  ! pour ecriture dans histxxx.nc
4028
4029  ! Get some variables from module fonte_neige_mod
4030  CALL fonte_neige_get_vars(pctsrf,  &
4031       zxfqcalving, zxfqfonte, zxffonte)
4032
4033
4034
4035
4036  !=============================================================
4037  ! Separation entre thermiques et non thermiques dans les sorties
4038  ! de fisrtilp
4039  !=============================================================
4040
4041  if (iflag_thermals>=1) then
4042     d_t_lscth=0.
4043     d_t_lscst=0.
4044     d_q_lscth=0.
4045     d_q_lscst=0.
4046     do k=1,klev
4047        do i=1,klon
4048           if (ptconvth(i,k)) then
4049              d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4050              d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4051           else
4052              d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4053              d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4054           endif
4055        enddo
4056     enddo
4057
4058     do i=1,klon
4059        plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4060        plul_th(i)=prfl(i,1)+psfl(i,1)
4061     enddo
4062  endif
4063
4064
4065  !On effectue les sorties:
4066
4067  CALL phys_output_write(itap, pdtphys, paprs, pphis,               &
4068       pplay, lmax_th, aerosol_couple,                 &
4069       ok_ade, ok_aie, ivap, new_aod, ok_sync,         &
4070       ptconv, read_climoz, clevSTD,                   &
4071       ptconvth, d_t, qx, d_qx, zmasse,                &
4072       flag_aerosol, flag_aerosol_strat, ok_cdnc)
4073
4074
4075
4076
4077  include "write_histday_seri.h"
4078
4079  include "write_paramLMDZ_phy.h"
4080
4081#endif
4082
4083  ! 22.03.04 END
4084  !
4085  !====================================================================
4086  ! Si c'est la fin, il faut conserver l'etat de redemarrage
4087  !====================================================================
4088  !
4089
4090  !        -----------------------------------------------------------------
4091  !        WSTATS: Saving statistics
4092  !        -----------------------------------------------------------------
4093  !        ("stats" stores and accumulates 8 key variables in file "stats.nc"
4094  !        which can later be used to make the statistic files of the run:
4095  !        "stats")          only possible in 3D runs !
4096
4097
4098  IF (callstats) THEN
4099
4100     call wstats(klon,o_psol%name,"Surface pressure","Pa" &
4101          ,2,paprs(:,1))
4102     call wstats(klon,o_tsol%name,"Surface temperature","K", &
4103          2,zxtsol)
4104     zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
4105     call wstats(klon,o_precip%name,"Precip Totale liq+sol", &
4106          "kg/(s*m2)",2,zx_tmp_fi2d)
4107     zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
4108     call wstats(klon,o_plul%name,"Large-scale Precip", &
4109          "kg/(s*m2)",2,zx_tmp_fi2d)
4110     zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
4111     call wstats(klon,o_pluc%name,"Convective Precip", &
4112          "kg/(s*m2)",2,zx_tmp_fi2d)
4113     call wstats(klon,o_sols%name,"Solar rad. at surf.", &
4114          "W/m2",2,solsw)
4115     call wstats(klon,o_soll%name,"IR rad. at surf.", &
4116          "W/m2",2,sollw)
4117     zx_tmp_fi2d(:) = topsw(:)-toplw(:)
4118     call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA", &
4119          "W/m2",2,zx_tmp_fi2d)
4120
4121
4122
4123     call wstats(klon,o_temp%name,"Air temperature","K", &
4124          3,t_seri)
4125     call wstats(klon,o_vitu%name,"Zonal wind","m.s-1", &
4126          3,u_seri)
4127     call wstats(klon,o_vitv%name,"Meridional wind", &
4128          "m.s-1",3,v_seri)
4129     call wstats(klon,o_vitw%name,"Vertical wind", &
4130          "m.s-1",3,omega)
4131     call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg", &
4132          3,q_seri)
4133
4134
4135
4136     IF(lafin) THEN
4137        write (*,*) "Writing stats..."
4138        call mkstats(ierr)
4139     ENDIF
4140
4141  ENDIF !if callstats
4142
4143  IF (lafin) THEN
4144     itau_phy = itau_phy + itap
4145     CALL phyredem ("restartphy.nc")
4146     !         open(97,form="unformatted",file="finbin")
4147     !         write(97) u_seri,v_seri,t_seri,q_seri
4148     !         close(97)
4149     !$OMP MASTER
4150     if (read_climoz >= 1) then
4151        if (is_mpi_root) then
4152           call nf95_close(ncid_climoz)
4153        end if
4154        deallocate(press_climoz) ! pointer
4155     end if
4156     !$OMP END MASTER
4157  ENDIF
4158
4159  !      first=.false.
4160
4161  RETURN
4162END SUBROUTINE physiq
4163FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
4164  IMPLICIT none
4165  !
4166  ! Calculer et imprimer l'eau totale. A utiliser pour verifier
4167  ! la conservation de l'eau
4168  !
4169  include "YOMCST.h"
4170  INTEGER klon,klev
4171  REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
4172  REAL aire(klon)
4173  REAL qtotal, zx, qcheck
4174  INTEGER i, k
4175  !
4176  zx = 0.0
4177  DO i = 1, klon
4178     zx = zx + aire(i)
4179  ENDDO
4180  qtotal = 0.0
4181  DO k = 1, klev
4182     DO i = 1, klon
4183        qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i) &
4184             *(paprs(i,k)-paprs(i,k+1))/RG
4185     ENDDO
4186  ENDDO
4187  !
4188  qcheck = qtotal/zx
4189  !
4190  RETURN
4191END FUNCTION qcheck
4192SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
4193  IMPLICIT none
4194  !
4195  ! Tranformer une variable de la grille physique a
4196  ! la grille d'ecriture
4197  !
4198  INTEGER nfield,nlon,iim,jjmp1, jjm
4199  REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
4200  !
4201  INTEGER i, n, ig
4202  !
4203  jjm = jjmp1 - 1
4204  DO n = 1, nfield
4205     DO i=1,iim
4206        ecrit(i,n) = fi(1,n)
4207        ecrit(i+jjm*iim,n) = fi(nlon,n)
4208     ENDDO
4209     DO ig = 1, nlon - 2
4210        ecrit(iim+ig,n) = fi(1+ig,n)
4211     ENDDO
4212  ENDDO
4213  RETURN
4214END SUBROUTINE gr_fi_ecrit
Note: See TracBrowser for help on using the repository browser.