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

Last change on this file since 2284 was 2278, checked in by fhourdin, 9 years ago

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