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

Last change on this file since 2171 was 2171, checked in by acozic, 10 years ago

There are some commits that we must not do just before holiday .... so be back to rev 2168

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