source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F @ 5017

Last change on this file since 5017 was 1399, checked in by idelkadi, 14 years ago

Modifications relatives a l'inclusion des nouveaux nuages d'Arnaud Jam

  • calltherm.F90 : remontee de variables suplementaires.

Passage de ztv, zpspsk, zthla, zthl en argument

  • fisrtilp.F : le programme de nuages avec appel a cloudth

Ajout des arguments : zqta, fraca, ztv, zpspsk, ztla, ztlh,
iflag_cldcon

  • physiq.F : appel modifie a calltherm et fisrtilp
  • thermcell_main.F90 : changement dans les flag_thermals_ed
  • thermcell_plume.F90 : version modifiee entrainement et detrainement

(on retrouve la precendente pour iflag_thermals_ed=10)

  • cloudth.F90 : PDF d'eau bi-gaussiennes
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 123.6 KB
Line 
1! $Id: physiq.F 1399 2010-06-08 08:29:11Z dcugnet $
2!
3c#define IO_DEBUG
4
5      SUBROUTINE physiq (nlon,nlev,
6     .            debut,lafin,jD_cur, jH_cur,pdtphys,
7     .            paprs,pplay,pphi,pphis,presnivs,clesphy0,
8     .            u,v,t,qx,
9     .            flxmass_w,
10     .            d_u, d_v, d_t, d_qx, d_ps
11     .            , dudyn
12     .            , PVteta)
13
14      USE ioipsl, only: histbeg, histvert, histdef, histend, histsync,
15     $     histwrite, ju2ymds, ymds2ju, ioget_year_len
16      USE comgeomphy
17      USE phys_cal_mod
18      USE write_field_phy
19      USE dimphy
20      USE infotrac
21      USE mod_grid_phy_lmdz
22      USE mod_phys_lmdz_para
23      USE iophy
24      USE misc_mod, mydebug=>debug
25      USE vampir
26      USE pbl_surface_mod, ONLY : pbl_surface
27      USE change_srf_frac_mod
28      USE surface_data,     ONLY : type_ocean, ok_veget
29      USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
30      USE phys_state_var_mod ! Variables sauvegardees de la physique
31      USE phys_output_var_mod ! Variables pour les ecritures des sorties
32      USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
33      USE phys_output_mod
34      use open_climoz_m, only: open_climoz ! ozone climatology from a file
35      use regr_pr_av_m, only: regr_pr_av
36      use netcdf95, only: nf95_close
37cIM for NMC files
38      use netcdf, only: nf90_fill_real
39      use mod_phys_lmdz_mpi_data, only: is_mpi_root
40      USE aero_mod
41      use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
42      use conf_phys_m, only: conf_phys
43      use radlwsw_m, only: radlwsw
44      USE control_mod
45
46
47      IMPLICIT none
48c======================================================================
49c
50c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
51c
52c Objet: Moniteur general de la physique du modele
53cAA      Modifications quant aux traceurs :
54cAA                  -  uniformisation des parametrisations ds phytrac
55cAA                  -  stockage des moyennes des champs necessaires
56cAA                     en mode traceur off-line
57c======================================================================
58c   CLEFS CPP POUR LES IO
59c   =====================
60#define histNMC
61c#define histISCCP
62c======================================================================
63c    modif   ( P. Le Van ,  12/10/98 )
64c
65c  Arguments:
66c
67c nlon----input-I-nombre de points horizontaux
68c nlev----input-I-nombre de couches verticales, doit etre egale a klev
69c debut---input-L-variable logique indiquant le premier passage
70c lafin---input-L-variable logique indiquant le dernier passage
71c jD_cur       -R-jour courant a l'appel de la physique (jour julien)
72c jH_cur       -R-heure courante a l'appel de la physique (jour julien)
73c pdtphys-input-R-pas d'integration pour la physique (seconde)
74c paprs---input-R-pression pour chaque inter-couche (en Pa)
75c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
76c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
77c pphis---input-R-geopotentiel du sol
78c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
79c u-------input-R-vitesse dans la direction X (de O a E) en m/s
80c v-------input-R-vitesse Y (de S a N) en m/s
81c t-------input-R-temperature (K)
82c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
83c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
84c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
85c flxmass_w -input-R- flux de masse verticale
86c d_u-----output-R-tendance physique de "u" (m/s/s)
87c d_v-----output-R-tendance physique de "v" (m/s/s)
88c d_t-----output-R-tendance physique de "t" (K/s)
89c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
90c d_ps----output-R-tendance physique de la pression au sol
91cIM
92c PVteta--output-R-vorticite potentielle a des thetas constantes
93c======================================================================
94#include "dimensions.h"
95      integer jjmp1
96      parameter (jjmp1=jjm+1-1/jjm)
97      integer iip1
98      parameter (iip1=iim+1)
99
100#include "regdim.h"
101#include "indicesol.h"
102#include "dimsoil.h"
103#include "clesphys.h"
104#include "temps.h"
105#include "iniprint.h"
106#include "thermcell.h"
107c======================================================================
108      LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE
109      PARAMETER (ok_cvl=.TRUE.)
110      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
111      PARAMETER (ok_gust=.FALSE.)
112      integer iflag_radia     ! active ou non le rayonnement (MPL)
113      save iflag_radia
114c$OMP THREADPRIVATE(iflag_radia)
115c======================================================================
116      LOGICAL check ! Verifier la conservation du modele en eau
117      PARAMETER (check=.FALSE.)
118      LOGICAL ok_stratus ! Ajouter artificiellement les stratus
119      PARAMETER (ok_stratus=.FALSE.)
120c======================================================================
121      REAL amn, amx
122      INTEGER igout
123c======================================================================
124c Clef controlant l'activation du cycle diurne:
125ccc      LOGICAL cycle_diurne
126ccc      PARAMETER (cycle_diurne=.FALSE.)
127c======================================================================
128c Modele thermique du sol, a activer pour le cycle diurne:
129ccc      LOGICAL soil_model
130ccc      PARAMETER (soil_model=.FALSE.)
131c======================================================================
132c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
133c le calcul du rayonnement est celle apres la precipitation des nuages.
134c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
135c la condensation et la precipitation. Cette cle augmente les impacts
136c radiatifs des nuages.
137ccc      LOGICAL new_oliq
138ccc      PARAMETER (new_oliq=.FALSE.)
139c======================================================================
140c Clefs controlant deux parametrisations de l'orographie:
141cc      LOGICAL ok_orodr
142ccc      PARAMETER (ok_orodr=.FALSE.)
143ccc      LOGICAL ok_orolf
144ccc      PARAMETER (ok_orolf=.FALSE.)
145c======================================================================
146      LOGICAL ok_journe ! sortir le fichier journalier
147      save ok_journe
148c$OMP THREADPRIVATE(ok_journe)
149c
150      LOGICAL ok_mensuel ! sortir le fichier mensuel
151      save ok_mensuel
152c$OMP THREADPRIVATE(ok_mensuel)
153c
154      LOGICAL ok_instan ! sortir le fichier instantane
155      save ok_instan
156c$OMP THREADPRIVATE(ok_instan)
157c
158      LOGICAL ok_LES ! sortir le fichier LES
159      save ok_LES                           
160c$OMP THREADPRIVATE(ok_LES)                 
161c
162      LOGICAL ok_region ! sortir le fichier regional
163      PARAMETER (ok_region=.FALSE.)
164c======================================================================
165      real weak_inversion(klon),dthmin(klon)
166      real seuil_inversion
167      save seuil_inversion
168c$OMP THREADPRIVATE(seuil_inversion)
169      integer iflag_ratqs
170      save iflag_ratqs
171c$OMP THREADPRIVATE(iflag_ratqs)
172      REAL lambda_th(klon,klev),zz,znum,zden
173      REAL wmax_th(klon)
174      REAL zmax_th(klon)
175      REAL tau_overturning_th(klon)
176
177      integer lmax_th(klon)
178      integer limbas(klon)
179      real ratqscth(klon,klev)
180      real ratqsdiff(klon,klev)
181      real zqsatth(klon,klev)
182
183c======================================================================
184c
185      INTEGER ivap          ! indice de traceurs pour vapeur d'eau
186      PARAMETER (ivap=1)
187      INTEGER iliq          ! indice de traceurs pour eau liquide
188      PARAMETER (iliq=2)
189
190c
191c
192c Variables argument:
193c
194      INTEGER nlon
195      INTEGER nlev
196      REAL, intent(in):: jD_cur, jH_cur
197
198      REAL pdtphys
199      LOGICAL debut, lafin
200      REAL paprs(klon,klev+1)
201      REAL pplay(klon,klev)
202      REAL pphi(klon,klev)
203      REAL pphis(klon)
204      REAL presnivs(klev)
205      REAL znivsig(klev)
206      real pir
207
208      REAL u(klon,klev)
209      REAL v(klon,klev)
210      REAL t(klon,klev),theta(klon,klev)
211      REAL qx(klon,klev,nqtot)
212      REAL flxmass_w(klon,klev)
213      REAL omega(klon,klev) ! vitesse verticale en Pa/s
214      REAL d_u(klon,klev)
215      REAL d_v(klon,klev)
216      REAL d_t(klon,klev)
217      REAL d_qx(klon,klev,nqtot)
218      REAL d_ps(klon)
219      real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
220!IM definition dynamique o_trac dans phys_output_open
221!      type(ctrl_out) :: o_trac(nqtot)
222c
223cIM Amip2 PV a theta constante
224c
225      INTEGER nbteta
226      PARAMETER(nbteta=3)
227      CHARACTER*3 ctetaSTD(nbteta)
228      DATA ctetaSTD/'350','380','405'/
229      SAVE ctetaSTD
230c$OMP THREADPRIVATE(ctetaSTD)
231      REAL rtetaSTD(nbteta)
232      DATA rtetaSTD/350., 380., 405./
233      SAVE rtetaSTD
234c$OMP THREADPRIVATE(rtetaSTD)     
235c
236      REAL PVteta(klon,nbteta)
237      REAL zx_tmp_3dte(iim,jjmp1,nbteta)
238c
239cMI Amip2 PV a theta constante
240
241cym      INTEGER klevp1, klevm1
242cym      PARAMETER(klevp1=klev+1,klevm1=klev-1)
243cym#include "raddim.h"
244c
245c
246cIM Amip2
247c variables a une pression donnee
248c
249      real rlevSTD(nlevSTD)
250      DATA rlevSTD/100000., 92500., 85000., 70000.,
251     .60000., 50000., 40000., 30000., 25000., 20000.,
252     .15000., 10000., 7000., 5000., 3000., 2000., 1000./
253      SAVE rlevstd
254c$OMP THREADPRIVATE(rlevstd)
255      CHARACTER*4 clevSTD(nlevSTD)
256      DATA clevSTD/'1000','925 ','850 ','700 ','600 ',
257     .'500 ','400 ','300 ','250 ','200 ','150 ','100 ',
258     .'70  ','50  ','30  ','20  ','10  '/
259      SAVE clevSTD
260c$OMP THREADPRIVATE(clevSTD)
261c
262      CHARACTER*4 bb2
263      CHARACTER*2 bb3
264c
265
266      real twriteSTD(klon,nlevSTD,nfiles)
267      real qwriteSTD(klon,nlevSTD,nfiles)
268      real rhwriteSTD(klon,nlevSTD,nfiles)
269      real phiwriteSTD(klon,nlevSTD,nfiles)
270      real uwriteSTD(klon,nlevSTD,nfiles)
271      real vwriteSTD(klon,nlevSTD,nfiles)
272      real wwriteSTD(klon,nlevSTD,nfiles)
273cIM for NMC files
274      REAL geo500(klon)
275      real :: rlevSTD3(nlevSTD3)
276      DATA rlevSTD3/85000., 50000., 25000./
277      SAVE rlevSTD3
278c$OMP THREADPRIVATE(rlevSTD3)
279      real :: rlevSTD8(nlevSTD8)
280      DATA rlevSTD8/100000., 85000., 70000., 50000., 25000., 10000.,
281     $     5000., 1000./
282      SAVE rlevSTD8
283c$OMP THREADPRIVATE(rlevSTD8)
284      real twriteSTD3(klon,nlevSTD3)
285      real qwriteSTD3(klon,nlevSTD3)
286      real rhwriteSTD3(klon,nlevSTD3)
287      real phiwriteSTD3(klon,nlevSTD3)
288      real uwriteSTD3(klon,nlevSTD3)
289      real vwriteSTD3(klon,nlevSTD3)
290      real wwriteSTD3(klon,nlevSTD3)
291c
292      real tnondefSTD8(klon,nlevSTD8)
293      real twriteSTD8(klon,nlevSTD8)
294      real qwriteSTD8(klon,nlevSTD8)
295      real rhwriteSTD8(klon,nlevSTD8)
296      real phiwriteSTD8(klon,nlevSTD8)
297      real uwriteSTD8(klon,nlevSTD8)
298      real vwriteSTD8(klon,nlevSTD8)
299      real wwriteSTD8(klon,nlevSTD8)
300c
301c plevSTD3 END
302c
303c nout : niveau de output des variables a une pression donnee
304      logical oknondef(klon,nlevSTD,nout)
305c
306c les produits uvSTD, vqSTD, .., T2STD sont calcules
307c a partir des valeurs instantannees toutes les 6 h
308c qui sont moyennees sur le mois
309c
310#include "radopt.h"
311c
312c
313c prw: precipitable water
314      real prw(klon)
315
316      REAL convliq(klon,klev)  ! eau liquide nuageuse convective
317      REAL convfra(klon,klev)  ! fraction nuageuse convective
318
319      REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
320      REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
321      REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
322      REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
323
324      INTEGER linv, kp1
325c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
326c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
327      REAL flwp(klon), fiwp(klon)
328      REAL flwc(klon,klev), fiwc(klon,klev)
329      REAL flwp_c(klon), fiwp_c(klon)
330      REAL flwc_c(klon,klev), fiwc_c(klon,klev)
331      REAL flwp_s(klon), fiwp_s(klon)
332      REAL flwc_s(klon,klev), fiwc_s(klon,klev)
333
334cIM ISCCP simulator v3.4
335c dans clesphys.h top_height, overlap
336cv3.4
337      INTEGER debug, debugcol
338cym      INTEGER npoints
339cym      PARAMETER(npoints=klon)
340c
341      INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
342      INTEGER nregISCtot
343      PARAMETER(nregISCtot=1)
344c
345c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
346c y compris pour 1 point
347c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
348c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
349      INTEGER imin_debut, nbpti
350      INTEGER jmin_debut, nbptj
351cIM parametres ISCCP BEG
352      INTEGER nbapp_isccp
353!     INTEGER nbapp_isccp,isccppas
354!     PARAMETER(isccppas=6) !appel du simulateurs tous les 6pas de temps de la physique
355!                           !i.e. toutes les 3 heures
356      INTEGER n
357      INTEGER ifreq_isccp(napisccp), freqin_pdt(napisccp)
358      DATA ifreq_isccp/3/
359      SAVE ifreq_isccp
360c$OMP THREADPRIVATE(ifreq_isccp)
361      CHARACTER*5 typinout(napisccp)
362      DATA typinout/'i3od'/
363      SAVE typinout
364c$OMP THREADPRIVATE(typinout)
365cIM verif boxptop BEG
366      CHARACTER*1 verticaxe(napisccp)
367      DATA verticaxe/'1'/
368      SAVE verticaxe
369c$OMP THREADPRIVATE(verticaxe)
370cIM verif boxptop END
371      INTEGER nvlev(napisccp)
372c     INTEGER nvlev
373      REAL t1, aa
374      REAL seed_re(klon,napisccp)
375cym !!!! A voir plus tard
376cym      INTEGER iphy(iim,jjmp1)
377cIM parametres ISCCP END
378c
379c ncol = nb. de sous-colonnes pour chaque maille du GCM
380c ncolmx = No. max. de sous-colonnes pour chaque maille du GCM
381c      INTEGER ncol(napisccp), ncolmx, seed(klon,napisccp)
382      INTEGER,SAVE :: ncol(napisccp)
383c$OMP THREADPRIVATE(ncol)
384      INTEGER ncolmx, seed(klon,napisccp)
385      REAL nbsunlit(nregISCtot,klon,napisccp)  !nbsunlit : moyenne de sunlit
386c     PARAMETER(ncolmx=1500)
387      PARAMETER(ncolmx=300)
388c
389cIM verif boxptop BEG
390      REAL vertlev(ncolmx,napisccp)
391cIM verif boxptop END
392c
393      REAL,SAVE :: tautab_omp(0:255),tautab(0:255)
394      INTEGER,SAVE :: invtau_omp(-20:45000),invtau(-20:45000)
395c$OMP THREADPRIVATE(tautab,invtau)
396      REAL emsfc_lw
397      PARAMETER(emsfc_lw=0.99)
398c     REAL    ran0                      ! type for random number fuction
399c
400      REAL cldtot(klon,klev)
401c variables de haut en bas pour le simulateur ISCCP
402      REAL dtau_s(klon,klev) !tau nuages startiformes
403      REAL dtau_c(klon,klev) !tau nuages convectifs
404      REAL dem_s(klon,klev)  !emissivite nuages startiformes
405      REAL dem_c(klon,klev)  !emissivite nuages convectifs
406c
407c variables de haut en bas pour le simulateur ISCCP
408      REAL pfull(klon,klev)
409      REAL phalf(klon,klev+1)
410      REAL qv(klon,klev)
411      REAL cc(klon,klev)
412      REAL conv(klon,klev)
413      REAL dtau_sH2B(klon,klev)
414      REAL dtau_cH2B(klon,klev)
415      REAL at(klon,klev)
416      REAL dem_sH2B(klon,klev)
417      REAL dem_cH2B(klon,klev)
418c
419      INTEGER kmax, lmax, lmax3
420      PARAMETER(kmax=8, lmax=8, lmax3=3)
421      INTEGER kmaxm1, lmaxm1
422      PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
423      INTEGER iimx7, jjmx7, jjmp1x7
424      PARAMETER(iimx7=iim*kmaxm1, jjmx7=jjm*lmaxm1,
425     .jjmp1x7=jjmp1*lmaxm1)
426c
427c output from ISCCP simulator
428      REAL fq_isccp(klon,kmaxm1,lmaxm1,napisccp)
429      REAL fq_is_true(klon,kmaxm1,lmaxm1,napisccp)
430      REAL totalcldarea(klon,napisccp)
431      REAL meanptop(klon,napisccp)
432      REAL meantaucld(klon,napisccp)
433      REAL boxtau(klon,ncolmx,napisccp)
434      REAL boxptop(klon,ncolmx,napisccp)
435      REAL zx_tmp_fi3d_bx(klon,ncolmx)
436      REAL zx_tmp_3d_bx(iim,jjmp1,ncolmx)
437c
438      REAL cld_fi3d(klon,lmax3)
439      REAL cld_3d(iim,jjmp1,lmax3)
440c
441      INTEGER iw, iwmax
442      REAL wmin, pas_w
443c     PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
444cIM 051005     PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
445      PARAMETER(wmin=-100.,pas_w=10.,iwmax=20)
446      REAL o500(klon)
447c
448
449c sorties ISCCP
450
451      integer nid_isccp
452      save nid_isccp       
453c$OMP THREADPRIVATE(nid_isccp)
454
455      REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
456      DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
457      SAVE zx_tau
458      DATA zx_pc/180., 310., 440., 560., 680., 800., 1000./
459      SAVE zx_pc
460c$OMP THREADPRIVATE(zx_tau,zx_pc)
461c cldtopres pression au sommet des nuages
462      REAL cldtopres(lmaxm1), cldtopres3(lmax3)
463      DATA cldtopres/180., 310., 440., 560., 680., 800., 1000./
464      DATA cldtopres3/440., 680., 1000./
465      SAVE cldtopres,cldtopres3
466c$OMP THREADPRIVATE(cldtopres,cldtopres3)
467cIM 051005 BEG
468      INTEGER komega, nhoriRD
469
470c taulev: numero du niveau de tau dans les sorties ISCCP
471      CHARACTER *4 taulev(kmaxm1)
472c     DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
473      DATA taulev/'tau0','tau1','tau2','tau3','tau4','tau5','tau6'/
474      CHARACTER *3 pclev(lmaxm1)
475      DATA pclev/'pc1','pc2','pc3','pc4','pc5','pc6','pc7'/
476      SAVE taulev,pclev
477c$OMP THREADPRIVATE(taulev,pclev)
478c
479c cnameisccp
480      CHARACTER *29 cnameisccp(lmaxm1,kmaxm1)
481cIM bad 151205     DATA cnameisccp/'pc< 50hPa, tau< 0.3',
482      DATA cnameisccp/'pc= 50-180hPa, tau< 0.3',
483     .                'pc= 180-310hPa, tau< 0.3',
484     .                'pc= 310-440hPa, tau< 0.3',
485     .                'pc= 440-560hPa, tau< 0.3',
486     .                'pc= 560-680hPa, tau< 0.3',
487     .                'pc= 680-800hPa, tau< 0.3',
488     .                'pc= 800-1000hPa, tau< 0.3',
489     .                'pc= 50-180hPa, tau= 0.3-1.3',
490     .                'pc= 180-310hPa, tau= 0.3-1.3',
491     .                'pc= 310-440hPa, tau= 0.3-1.3',
492     .                'pc= 440-560hPa, tau= 0.3-1.3',
493     .                'pc= 560-680hPa, tau= 0.3-1.3',
494     .                'pc= 680-800hPa, tau= 0.3-1.3',
495     .                'pc= 800-1000hPa, tau= 0.3-1.3',
496     .                'pc= 50-180hPa, tau= 1.3-3.6',
497     .                'pc= 180-310hPa, tau= 1.3-3.6',
498     .                'pc= 310-440hPa, tau= 1.3-3.6',
499     .                'pc= 440-560hPa, tau= 1.3-3.6',
500     .                'pc= 560-680hPa, tau= 1.3-3.6',
501     .                'pc= 680-800hPa, tau= 1.3-3.6',
502     .                'pc= 800-1000hPa, tau= 1.3-3.6',
503     .                'pc= 50-180hPa, tau= 3.6-9.4',
504     .                'pc= 180-310hPa, tau= 3.6-9.4',
505     .                'pc= 310-440hPa, tau= 3.6-9.4',
506     .                'pc= 440-560hPa, tau= 3.6-9.4',
507     .                'pc= 560-680hPa, tau= 3.6-9.4',
508     .                'pc= 680-800hPa, tau= 3.6-9.4',
509     .                'pc= 800-1000hPa, tau= 3.6-9.4',
510     .                'pc= 50-180hPa, tau= 9.4-23',
511     .                'pc= 180-310hPa, tau= 9.4-23',
512     .                'pc= 310-440hPa, tau= 9.4-23',
513     .                'pc= 440-560hPa, tau= 9.4-23',
514     .                'pc= 560-680hPa, tau= 9.4-23',
515     .                'pc= 680-800hPa, tau= 9.4-23',
516     .                'pc= 800-1000hPa, tau= 9.4-23',
517     .                'pc= 50-180hPa, tau= 23-60',
518     .                'pc= 180-310hPa, tau= 23-60',
519     .                'pc= 310-440hPa, tau= 23-60',
520     .                'pc= 440-560hPa, tau= 23-60',
521     .                'pc= 560-680hPa, tau= 23-60',
522     .                'pc= 680-800hPa, tau= 23-60',
523     .                'pc= 800-1000hPa, tau= 23-60',
524     .                'pc= 50-180hPa, tau> 60.',
525     .                'pc= 180-310hPa, tau> 60.',
526     .                'pc= 310-440hPa, tau> 60.',
527     .                'pc= 440-560hPa, tau> 60.',
528     .                'pc= 560-680hPa, tau> 60.',
529     .                'pc= 680-800hPa, tau> 60.',
530     .                'pc= 800-1000hPa, tau> 60.'/
531       SAVE cnameisccp
532c$OMP THREADPRIVATE(cnameisccp)
533c
534c     REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
535c     INTEGER nhorix7
536cIM: region='3d' <==> sorties en global
537      CHARACTER*3 region
538      PARAMETER(region='3d')
539c
540cIM ISCCP simulator v3.4
541c
542      logical ok_hf
543c
544      integer nid_hf, nid_hf3d
545      save ok_hf, nid_hf, nid_hf3d
546c$OMP THREADPRIVATE(ok_hf, nid_hf, nid_hf3d)
547c  QUESTION : noms de variables ?
548
549      INTEGER        longcles
550      PARAMETER    ( longcles = 20 )
551      REAL clesphy0( longcles      )
552c
553c Variables propres a la physique
554      INTEGER itap
555      SAVE itap                   ! compteur pour la physique
556c$OMP THREADPRIVATE(itap)
557c
558      real slp(klon) ! sea level pressure
559c
560      REAL fevap(klon,nbsrf)
561      REAL fluxlat(klon,nbsrf)
562c
563      REAL qsol(klon)
564      REAL,save ::  solarlong0
565c$OMP THREADPRIVATE(solarlong0)
566
567c
568c  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
569c
570cIM 141004     REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
571      REAL zulow(klon),zvlow(klon)
572c
573      INTEGER igwd,idx(klon),itest(klon)
574c
575      REAL agesno(klon,nbsrf)
576c
577c      REAL,allocatable,save :: run_off_lic_0(:)
578cc$OMP THREADPRIVATE(run_off_lic_0)
579cym      SAVE run_off_lic_0
580cKE43
581c Variables liees a la convection de K. Emanuel (sb):
582c
583      REAL bas, top             ! cloud base and top levels
584      SAVE bas
585      SAVE top
586c$OMP THREADPRIVATE(bas, top)
587
588      REAL wdn(klon), tdn(klon), qdn(klon)
589c
590c=================================================================================================
591cCR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
592c Variables liées à la poche froide (jyg)
593
594      REAL mip(klon,klev)  ! mass flux shed by the adiab ascent at each level
595      REAL Vprecip(klon,klev+1)   ! precipitation vertical profile
596c
597      REAL wape_prescr, fip_prescr
598      INTEGER it_wape_prescr
599      SAVE wape_prescr, fip_prescr, it_wape_prescr
600c$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
601c
602c variables supplementaires de concvl
603      REAL Tconv(klon,klev)
604      REAL ment(klon,klev,klev),sij(klon,klev,klev)
605      REAL dd_t(klon,klev),dd_q(klon,klev)
606
607      real, save :: alp_bl_prescr=0.
608      real, save :: ale_bl_prescr=0.
609
610      real, save :: ale_max=100.
611      real, save :: alp_max=2.
612
613c$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
614c$OMP THREADPRIVATE(ale_max,alp_max)
615
616      real ale_wake(klon)
617      real alp_wake(klon)
618cRC
619c Variables liées à la poche froide (jyg et rr)
620c Version diagnostique pour l'instant : pas de rétroaction sur la convection
621
622      REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
623
624      REAL wake_dth(klon,klev)        ! wake : temp pot difference
625
626      REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to Gravity Wave (/s)
627      REAL wake_omgbdth(klon,klev)    ! Wake : flux of Delta_Theta transported by LS omega
628      REAL wake_dp_omgb(klon,klev)    ! Wake : vertical gradient of large scale omega
629      REAL wake_dtKE(klon,klev)       ! Wake : differential heating (wake - unpertubed) CONV
630      REAL wake_dqKE(klon,klev)       ! Wake : differential moistening (wake - unpertubed) CONV
631      REAL wake_dtPBL(klon,klev)      ! Wake : differential heating (wake - unpertubed) PBL
632      REAL wake_dqPBL(klon,klev)      ! Wake : differential moistening (wake - unpertubed) PBL
633      REAL wake_omg(klon,klev)        ! Wake : velocity difference (wake - unpertubed)
634      REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
635      REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
636      REAL wake_spread(klon,klev)     ! spreading term in wake_delt
637c
638cpourquoi y'a pas de save??
639      REAL wake_h(klon)               ! Wake : hauteur de la poche froide
640c
641      INTEGER wake_k(klon)            ! Wake sommet
642c
643      REAL t_undi(klon,klev)               ! temperature moyenne dans la zone non perturbee
644      REAL q_undi(klon,klev)               ! humidite moyenne dans la zone non perturbee
645c
646cjyg
647ccc      REAL wake_pe(klon)              ! Wake potential energy - WAPE
648
649      REAL wake_gfl(klon)             ! Gust Front Length
650      REAL wake_dens(klon)
651c
652c
653      REAL dt_dwn(klon,klev)
654      REAL dq_dwn(klon,klev)
655      REAL wdt_PBL(klon,klev)
656      REAL udt_PBL(klon,klev)
657      REAL wdq_PBL(klon,klev)
658      REAL udq_PBL(klon,klev)
659      REAL M_dwn(klon,klev)
660      REAL M_up(klon,klev)
661      REAL dt_a(klon,klev)
662      REAL dq_a(klon,klev)
663      REAL, SAVE :: alp_offset
664c$OMP THREADPRIVATE(alp_offset)
665
666c
667cRR:fin declarations poches froides
668c=======================================================================================================
669
670      REAL zw2(klon,klev+1)
671      REAL fraca(klon,klev+1)       
672      REAL ztv(klon,klev)
673      REAL zpspsk(klon,klev)
674      REAL ztla(klon,klev)
675      REAL zthl(klon,klev)
676
677c Variables locales pour la couche limite (al1):
678c
679cAl1      REAL pblh(klon)           ! Hauteur de couche limite
680cAl1      SAVE pblh
681c34EK
682c
683c Variables locales:
684c
685      REAL cdragh(klon) ! drag coefficient pour T and Q
686      REAL cdragm(klon) ! drag coefficient pour vent
687cAA
688cAA  Pour phytrac
689cAA
690      REAL coefh(klon,klev)     ! coef d'echange pour phytrac, valable pour 2<=k<=klev
691      REAL u1(klon)             ! vents dans la premiere couche U
692      REAL v1(klon)             ! vents dans la premiere couche V
693
694      REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
695
696c@$$      LOGICAL offline           ! Controle du stockage ds "physique"
697c@$$      PARAMETER (offline=.false.)
698c@$$      INTEGER physid
699      REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
700      REAL frac_nucl(klon,klev) ! idem (nucleation)
701      INTEGER       :: iii
702      REAL          :: calday
703
704cIM cf FH pour Tiedtke 080604
705      REAL rain_tiedtke(klon),snow_tiedtke(klon)
706c
707cIM 050204 END
708      REAL evap(klon), devap(klon) ! evaporation et sa derivee
709      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
710
711      REAL bils(klon) ! bilan de chaleur au sol
712      REAL wfbilo(klon,nbsrf) ! bilan d'eau, pour chaque
713C                             ! type de sous-surface et pondere par la fraction
714      REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
715C                             ! type de sous-surface et pondere par la fraction
716      REAL slab_wfbils(klon)  ! bilan de chaleur au sol pour le cas de slab, sur les points d'ocean
717
718      REAL fder(klon)         
719      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
720      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
721      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
722      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
723c
724      REAL frugs(klon,nbsrf)
725      REAL zxrugs(klon) ! longueur de rugosite
726c
727c Conditions aux limites
728c
729!
730      REAL :: day_since_equinox
731! Date de l'equinoxe de printemps
732      INTEGER, parameter :: mth_eq=3, day_eq=21
733      REAL :: jD_eq
734
735      LOGICAL, parameter :: new_orbit = .true.
736
737c
738      INTEGER lmt_pas
739      SAVE lmt_pas                ! frequence de mise a jour
740c$OMP THREADPRIVATE(lmt_pas)
741      real zmasse(klon, llm)
742C     (column-density of mass of air in a cell, in kg m-2)
743      real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
744
745cIM sorties
746      REAL un_jour
747      PARAMETER(un_jour=86400.)
748c======================================================================
749c
750c Declaration des procedures appelees
751c
752      EXTERNAL angle     ! calculer angle zenithal du soleil
753      EXTERNAL alboc     ! calculer l'albedo sur ocean
754      EXTERNAL ajsec     ! ajustement sec
755      EXTERNAL conlmd    ! convection (schema LMD)
756cKE43
757      EXTERNAL conema3  ! convect4.3
758      EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
759cAA
760      EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
761c                          ! stockage des coefficients necessaires au
762c                          ! lessivage OFF-LINE et ON-LINE
763      EXTERNAL hgardfou  ! verifier les temperatures
764      EXTERNAL nuage     ! calculer les proprietes radiatives
765CC      EXTERNAL o3cm      ! initialiser l'ozone
766      EXTERNAL orbite    ! calculer l'orbite terrestre
767      EXTERNAL phyetat0  ! lire l'etat initial de la physique
768      EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
769      EXTERNAL suphel    ! initialiser certaines constantes
770      EXTERNAL transp    ! transport total de l'eau et de l'energie
771      EXTERNAL ecribina  ! ecrire le fichier binaire global
772      EXTERNAL ecribins  ! ecrire le fichier binaire global
773      EXTERNAL ecrirega  ! ecrire le fichier binaire regional
774      EXTERNAL ecriregs  ! ecrire le fichier binaire regional
775cIM
776      EXTERNAL haut2bas  !variables de haut en bas
777      INTEGER lnblnk1
778      EXTERNAL lnblnk1   !enleve les blancs a la fin d'une variable de type
779                         !caracter
780      EXTERNAL ini_undefSTD  !initialise a 0 une variable a 1 niveau de pression
781      EXTERNAL undefSTD      !somme les valeurs definies d'1 var a 1 niveau de pression
782c     EXTERNAL moy_undefSTD  !moyenne d'1 var a 1 niveau de pression
783c     EXTERNAL moyglo_aire   !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
784c                            !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
785c
786c Variables locales
787c
788      REAL rhcl(klon,klev)    ! humiditi relative ciel clair
789      REAL dialiq(klon,klev)  ! eau liquide nuageuse
790      REAL diafra(klon,klev)  ! fraction nuageuse
791      REAL cldliq(klon,klev)  ! eau liquide nuageuse
792      REAL cldfra(klon,klev)  ! fraction nuageuse
793      REAL cldtau(klon,klev)  ! epaisseur optique
794      REAL cldemi(klon,klev)  ! emissivite infrarouge
795c
796CXXX PB
797      REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
798      REAL fluxt(klon,klev, nbsrf)   ! flux turbulent de chaleur
799      REAL fluxu(klon,klev, nbsrf)   ! flux turbulent de vitesse u
800      REAL fluxv(klon,klev, nbsrf)   ! flux turbulent de vitesse v
801c
802      REAL zxfluxt(klon, klev)
803      REAL zxfluxq(klon, klev)
804      REAL zxfluxu(klon, klev)
805      REAL zxfluxv(klon, klev)
806CXXX
807c
808      REAL fsollw(klon, nbsrf)   ! bilan flux IR pour chaque sous surface
809      REAL fsolsw(klon, nbsrf)   ! flux solaire absorb. pour chaque sous surface
810c Le rayonnement n'est pas calcule tous les pas, il faut donc
811c                      sauvegarder les sorties du rayonnement
812cym      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
813cym      SAVE  sollwdownclr, toplwdown, toplwdownclr
814cym      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
815c
816      INTEGER itaprad
817      SAVE itaprad
818c$OMP THREADPRIVATE(itaprad)
819c
820      REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
821      REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
822c
823      REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
824      REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
825c
826      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
827      REAL zxsnow_dummy(klon)
828c
829      REAL dist, rmu0(klon), fract(klon)
830      REAL zdtime, zlongi
831c
832      CHARACTER*2 str2
833      CHARACTER*2 iqn
834c
835      REAL qcheck
836      REAL z_avant(klon), z_apres(klon), z_factor(klon)
837      LOGICAL zx_ajustq
838c
839      REAL za, zb
840      REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
841      real zqsat(klon,klev)
842      INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq, iff
843      REAL t_coup
844      PARAMETER (t_coup=234.0)
845c
846      REAL zphi(klon,klev)
847cym A voir plus tard !!
848cym      REAL zx_relief(iim,jjmp1)
849cym      REAL zx_aire(iim,jjmp1)
850c
851c Grandeurs de sorties
852      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
853      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
854      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
855      REAL s_trmb3(klon)
856cKE43
857c Variables locales pour la convection de K. Emanuel (sb):
858c
859      REAL upwd(klon,klev)      ! saturated updraft mass flux
860      REAL dnwd(klon,klev)      ! saturated downdraft mass flux
861      REAL dnwd0(klon,klev)     ! unsaturated downdraft mass flux
862      REAL tvp(klon,klev)       ! virtual temp of lifted parcel
863      CHARACTER*40 capemaxcels  !max(CAPE)
864
865      REAL rflag(klon)          ! flag fonctionnement de convect
866      INTEGER iflagctrl(klon)          ! flag fonctionnement de convect
867c -- convect43:
868      INTEGER ntra              ! nb traceurs pour convect4.3
869      REAL pori_con(klon)    ! pressure at the origin level of lifted parcel
870      REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
871      REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
872      REAL dplcldt(klon), dplcldr(klon)
873c?     .     condm_con(klon,klev),conda_con(klon,klev),
874c?     .     mr_con(klon,klev),ep_con(klon,klev)
875c?     .    ,sadiab(klon,klev),wadiab(klon,klev)
876c --
877c34EK
878c
879c Variables du changement
880c
881c con: convection
882c lsc: condensation a grande echelle (Large-Scale-Condensation)
883c ajs: ajustement sec
884c eva: evaporation de l'eau liquide nuageuse
885c vdf: couche limite (Vertical DiFfusion)
886      REAL rneb(klon,klev)
887
888! tendance nulles
889      REAL du0(klon,klev),dv0(klon,klev),dq0(klon,klev),dql0(klon,klev)
890
891c
892*********************************************************
893*     declarations
894     
895*********************************************************
896cIM 081204 END
897c
898      REAL pmfu(klon,klev), pmfd(klon,klev)
899      REAL pen_u(klon,klev), pen_d(klon,klev)
900      REAL pde_u(klon,klev), pde_d(klon,klev)
901      INTEGER kcbot(klon), kctop(klon), kdtop(klon)
902      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
903      REAL prfl(klon,klev+1), psfl(klon,klev+1)
904c
905      REAL rain_lsc(klon)
906      REAL snow_lsc(klon)
907c
908      REAL ratqss(klon,klev),ratqsc(klon,klev)
909      real ratqsbas,ratqshaut,tau_ratqs
910      save ratqsbas,ratqshaut,tau_ratqs
911c$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
912      real zpt_conv(klon,klev)
913
914c Parametres lies au nouveau schema de nuages (SB, PDF)
915      real fact_cldcon
916      real facttemps
917      logical ok_newmicro
918      save ok_newmicro
919      real ref_liq(klon,klev), ref_ice(klon,klev)
920c$OMP THREADPRIVATE(ok_newmicro)
921      save fact_cldcon,facttemps
922c$OMP THREADPRIVATE(fact_cldcon,facttemps)
923      real facteur
924
925      integer iflag_cldcon
926      save iflag_cldcon
927c$OMP THREADPRIVATE(iflag_cldcon)
928      logical ptconv(klon,klev)
929cIM cf. AM 081204 BEG
930      logical ptconvth(klon,klev)
931cIM cf. AM 081204 END
932c
933c Variables liees a l'ecriture de la bande histoire physique
934c
935c======================================================================
936c
937cIM cf. AM 081204 BEG
938c   declarations pour sortir sur une sous-region
939      integer imin_ins,imax_ins,jmin_ins,jmax_ins
940      save imin_ins,imax_ins,jmin_ins,jmax_ins
941c$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
942c      real lonmin_ins,lonmax_ins,latmin_ins
943c     s  ,latmax_ins
944c     data lonmin_ins,lonmax_ins,latmin_ins
945c    s  ,latmax_ins/
946c valeurs initiales     s   -5.,20.,41.,55./   
947c    s   100.,130.,-20.,20./
948c    s   -180.,180.,-90.,90./
949c======================================================================
950cIM cf. AM 081204 END
951
952c
953      integer itau_w   ! pas de temps ecriture = itap + itau_phy
954c
955c
956c Variables locales pour effectuer les appels en serie
957c
958      REAL zx_rh(klon,klev)
959cIM RH a 2m (la surface)
960      REAL rh2m(klon), qsat2m(klon)
961      REAL tpot(klon), tpote(klon)
962      REAL Lheat
963
964      INTEGER        length
965      PARAMETER    ( length = 100 )
966      REAL tabcntr0( length       )
967c
968      INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
969cIM
970      INTEGER ndex2d1(iwmax)
971c
972cIM AMIP2 BEG
973      REAL moyglo, mountor
974cIM 141004 BEG
975      REAL zustrdr(klon), zvstrdr(klon)
976      REAL zustrli(klon), zvstrli(klon)
977      REAL zustrph(klon), zvstrph(klon)
978      REAL zustrhi(klon), zvstrhi(klon)
979      REAL aam, torsfc
980cIM 141004 END
981cIM 190504 BEG
982      INTEGER ij, imp1jmp1
983      PARAMETER(imp1jmp1=(iim+1)*jjmp1)
984cym A voir plus tard
985      REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
986      REAL padyn(iim+1,jjmp1,klev+1)
987      REAL dudyn(iim+1,jjmp1,klev)
988      REAL rlatdyn(iim+1,jjmp1)
989cIM 190504 END
990      LOGICAL ok_msk
991      REAL msk(klon)
992cIM
993      REAL airetot, pi
994cym A voir plus tard
995cym      REAL zm_wo(jjmp1, klev)
996cIM AMIP2 END
997c
998      REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
999      REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
1000c#ifdef histNMC
1001cym   A voir plus tard !!!!
1002cym      REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
1003      REAL zx_tmp_fiNC(klon,nlevSTD)
1004c#endif
1005      REAL(KIND=8) zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
1006      REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
1007      REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
1008cIM for NMC files
1009      REAL missing_val
1010      REAL, SAVE :: freq_moyNMC(nout)
1011c$OMP THREADPRIVATE(freq_moyNMC)
1012c
1013      INTEGER nid_day, nid_mth, nid_ins, nid_mthnmc, nid_daynmc
1014      INTEGER nid_hfnmc, nid_day_seri, nid_ctesGCM
1015      SAVE nid_day, nid_mth, nid_ins, nid_mthnmc, nid_daynmc
1016      SAVE nid_hfnmc, nid_day_seri, nid_ctesGCM
1017c$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins)
1018c$OMP THREADPRIVATE(nid_mthnmc, nid_daynmc, nid_hfnmc)
1019c$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
1020c
1021cIM 280405 BEG
1022      INTEGER nid_bilKPins, nid_bilKPave
1023      SAVE nid_bilKPins, nid_bilKPave
1024c$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
1025c
1026      REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
1027      REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
1028      REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
1029      REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
1030c
1031cIM 280405 END
1032c
1033      INTEGER nhori, nvert, nvert1, nvert3
1034      REAL zsto, zsto1, zsto2
1035      REAL zstophy, zstorad, zstohf, zstoday, zstomth, zout
1036      REAL zcals(napisccp), zcalh(napisccp), zoutj(napisccp)
1037      REAL zout_isccp(napisccp)
1038      SAVE zcals, zcalh, zoutj, zout_isccp
1039c$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
1040
1041      real zjulian
1042      save zjulian
1043c$OMP THREADPRIVATE(zjulian)
1044
1045      character*20 modname
1046      character*80 abort_message
1047      logical ok_sync
1048      real date0
1049      integer idayref
1050
1051C essai writephys
1052      integer fid_day, fid_mth, fid_ins
1053      parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
1054      integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
1055      parameter (prof2d_on = 1, prof3d_on = 2,
1056     .           prof2d_av = 3, prof3d_av = 4)
1057      character*30 nom_fichier
1058      character*10 varname
1059      character*40 vartitle
1060      character*20 varunits
1061C     Variables liees au bilan d'energie et d'enthalpi
1062      REAL ztsol(klon)
1063      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1064     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1065      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1066     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1067c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot,
1068c$OMP+              h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
1069      REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
1070      REAL      d_h_vcol_phy
1071      REAL      fs_bound, fq_bound
1072      SAVE      d_h_vcol_phy
1073c$OMP THREADPRIVATE(d_h_vcol_phy)
1074      REAL      zero_v(klon)
1075      CHARACTER*15 ztit
1076      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
1077      SAVE      ip_ebil
1078      DATA      ip_ebil/0/
1079c$OMP THREADPRIVATE(ip_ebil)
1080      INTEGER   if_ebil ! level for energy conserv. dignostics
1081      SAVE      if_ebil
1082c$OMP THREADPRIVATE(if_ebil)
1083c+jld ec_conser
1084      REAL ZRCPD
1085c-jld ec_conser
1086      REAL t2m(klon,nbsrf)  ! temperature a 2m
1087      REAL q2m(klon,nbsrf)  ! humidite a 2m
1088
1089cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
1090      REAL zt2m(klon), zq2m(klon)             !temp., hum. 2m moyenne s/ 1 maille
1091      REAL zu10m(klon), zv10m(klon)           !vents a 10m moyennes s/1 maille
1092      CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
1093      CHARACTER*40 tinst, tave, typeval
1094      REAL cldtaupi(klon,klev)  ! Cloud optical thickness for pre-industrial (pi) aerosols
1095
1096      REAL re(klon, klev)       ! Cloud droplet effective radius
1097      REAL fl(klon, klev)  ! denominator of re
1098
1099      REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
1100
1101      ! Aerosol optical properties
1102      CHARACTER*4, DIMENSION(naero_grp) :: rfname
1103      REAL, DIMENSION(klon)          :: aerindex     ! POLDER aerosol index
1104      REAL, DIMENSION(klon,klev)     :: mass_solu_aero    ! total mass concentration for all soluble aerosols[ug/m3]
1105      REAL, DIMENSION(klon,klev)     :: mass_solu_aero_pi ! - " - (pre-industrial value)
1106      INTEGER :: naero ! aerosol species
1107
1108      ! Parameters
1109      LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not
1110      REAL bl95_b0, bl95_b1   ! Parameter in Boucher and Lohmann (1995)
1111      SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
1112c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
1113      LOGICAL, SAVE :: aerosol_couple ! true  : calcul des aerosols dans INCA
1114                                      ! false : lecture des aerosol dans un fichier
1115c$OMP THREADPRIVATE(aerosol_couple)   
1116      INTEGER, SAVE :: flag_aerosol
1117c$OMP THREADPRIVATE(flag_aerosol)
1118      LOGICAL, SAVE :: new_aod
1119c$OMP THREADPRIVATE(new_aod)
1120   
1121c
1122c Declaration des constantes et des fonctions thermodynamiques
1123c
1124      LOGICAL,SAVE :: first=.true.
1125c$OMP THREADPRIVATE(first)
1126
1127      integer iunit
1128
1129      integer, save::  read_climoz ! read ozone climatology
1130C     Allowed values are 0, 1 and 2
1131C     0: do not read an ozone climatology
1132C     1: read a single ozone climatology that will be used day and night
1133C     2: read two ozone climatologies, the average day and night
1134C     climatology and the daylight climatology
1135
1136      integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
1137
1138      real, pointer, save:: press_climoz(:)
1139!     edges of pressure intervals for ozone climatologies, in Pa, in strictly
1140!     ascending order
1141
1142      integer, save:: co3i = 0
1143!     time index in NetCDF file of current ozone fields
1144c$OMP THREADPRIVATE(co3i)
1145
1146      integer ro3i
1147!     required time index in NetCDF file for the ozone fields, between 1
1148!     and 360
1149
1150#include "YOMCST.h"
1151#include "YOETHF.h"
1152#include "FCTTRE.h"
1153cIM 100106 BEG : pouvoir sortir les ctes de la physique
1154#include "conema3.h"
1155#include "fisrtilp.h"
1156#include "nuage.h"
1157#include "compbl.h"
1158cIM 100106 END : pouvoir sortir les ctes de la physique
1159c
1160!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1161c Declarations pour Simulateur COSP
1162c============================================================
1163      real :: mr_ozone(klon,klev)
1164cIM for NMC files
1165      missing_val=nf90_fill_real
1166c======================================================================
1167! Gestion calendrier : mise a jour du module phys_cal_mod
1168!
1169      CALL phys_cal_update(jD_cur,jH_cur)
1170
1171c======================================================================
1172! Ecriture eventuelle d'un profil verticale en entree de la physique.
1173! Utilise notamment en 1D mais peut etre active egalement en 3D
1174! en imposant la valeur de igout.
1175c======================================================================
1176
1177      if (prt_level.ge.1) then
1178          igout=klon/2+1/klon
1179         write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1180         write(lunout,*)
1181     s 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1182         write(lunout,*)
1183     s  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1184
1185         write(lunout,*) 'paprs, play, phi, u, v, t'
1186         do k=1,klev
1187            write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
1188     s   u(igout,k),v(igout,k),t(igout,k)
1189         enddo
1190         write(lunout,*) 'ovap (g/kg),  oliq (g/kg)'
1191         do k=1,klev
1192            write(lunout,*) qx(igout,k,1)*1000,qx(igout,k,2)*1000.
1193         enddo
1194      endif
1195
1196c======================================================================
1197
1198cym => necessaire pour iflag_con != 2   
1199      pmfd(:,:) = 0.
1200      pen_u(:,:) = 0.
1201      pen_d(:,:) = 0.
1202      pde_d(:,:) = 0.
1203      pde_u(:,:) = 0.
1204      aam=0.
1205
1206      torsfc=0.
1207      forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1208
1209      if (first) then
1210     
1211cCR:nvelles variables convection/poches froides
1212     
1213      print*, '================================================='
1214      print*, 'Allocation des variables locales et sauvegardees'
1215      call phys_local_var_init
1216c
1217      pasphys=pdtphys
1218c     appel a la lecture du run.def physique
1219      call conf_phys(ok_journe, ok_mensuel,
1220     .     ok_instan, ok_hf,
1221     .     ok_LES,
1222     .     solarlong0,seuil_inversion,
1223     .     fact_cldcon, facttemps,ok_newmicro,iflag_radia,
1224     .     iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,
1225     .     ok_ade, ok_aie, aerosol_couple,
1226     .     flag_aerosol, new_aod,
1227     .     bl95_b0, bl95_b1,
1228     .     iflag_thermals,nsplit_thermals,tau_thermals,
1229     .     iflag_thermals_ed,iflag_thermals_optflux,
1230c     nv flags pour la convection et les poches froides
1231     .     iflag_coupl,iflag_clos,iflag_wake, read_climoz,
1232     &     alp_offset)
1233      call phys_state_var_init(read_climoz)
1234      call phys_output_var_init
1235
1236      print*, '================================================='
1237cIM for NMC files
1238cIM freq_moyNMC = frequences auxquelles on moyenne les champs accumules
1239cIM               sur les niveaux de pression standard du NMC
1240      DO n=1, nout
1241       freq_moyNMC(n)=freq_outNMC(n)/freq_calNMC(n)
1242      ENDDO
1243c
1244cIM beg
1245          dnwd0=0.0
1246          ftd=0.0
1247          fqd=0.0
1248          cin=0.
1249cym Attention pbase pas initialise dans concvl !!!!
1250          pbase=0
1251          paire_ter(:)=0.   
1252cIM 180608
1253c         pmflxr=0.
1254c         pmflxs=0.
1255        itau_con=0
1256        first=.false.
1257
1258      endif  ! first
1259
1260       modname = 'physiq'
1261cIM
1262      IF (ip_ebil_phy.ge.1) THEN
1263        DO i=1,klon
1264          zero_v(i)=0.
1265        END DO
1266      END IF
1267      ok_sync=.TRUE.
1268
1269      IF (debut) THEN
1270         CALL suphel ! initialiser constantes et parametres phys.
1271      ENDIF
1272
1273      if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1274
1275
1276c======================================================================
1277! Gestion calendrier : mise a jour du module phys_cal_mod
1278!
1279c     CALL phys_cal_update(jD_cur,jH_cur)
1280
1281c
1282c Si c'est le debut, il faut initialiser plusieurs choses
1283c          ********
1284c
1285       IF (debut) THEN
1286!rv
1287cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1288cde la convection a partir des caracteristiques du thermique
1289         wght_th(:,:)=1.
1290         lalim_conv(:)=1
1291cRC
1292         u10m(:,:)=0.
1293         v10m(:,:)=0.
1294         rain_con(:)=0.
1295         snow_con(:)=0.
1296         topswai(:)=0.
1297         topswad(:)=0.
1298         solswai(:)=0.
1299         solswad(:)=0.
1300
1301         lambda_th(:,:)=0.
1302         wmax_th(:)=0.
1303         tau_overturning_th(:)=0.
1304
1305         IF (config_inca /= 'none') THEN
1306            ! jg : initialisation jusqu'au ces variables sont dans restart
1307            ccm(:,:,:) = 0.
1308            tau_aero(:,:,:,:) = 0.
1309            piz_aero(:,:,:,:) = 0.
1310            cg_aero(:,:,:,:) = 0.
1311         END IF
1312
1313         rnebcon0(:,:) = 0.0
1314         clwcon0(:,:) = 0.0
1315         rnebcon(:,:) = 0.0
1316         clwcon(:,:) = 0.0
1317
1318cIM     
1319         IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1320c
1321      print*,'iflag_coupl,iflag_clos,iflag_wake',
1322     .   iflag_coupl,iflag_clos,iflag_wake
1323      print*,'CYCLE_DIURNE', cycle_diurne
1324c
1325      IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
1326         abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
1327         CALL abort_gcm (modname,abort_message,1)
1328      ENDIF
1329c
1330      IF(ok_isccp.AND.iflag_con.LE.2) THEN
1331         abort_message = 'ISCCP-like outputs may be available for KE
1332     .(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
1333         CALL abort_gcm (modname,abort_message,1)
1334      ENDIF
1335c
1336c Initialiser les compteurs:
1337c
1338         itap    = 0
1339         itaprad = 0
1340
1341!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1342!! Un petit travail à faire ici.
1343!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1344
1345         if (iflag_pbl>1) then
1346             PRINT*, "Using method MELLOR&YAMADA"
1347         endif
1348
1349!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1350! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1351! dyn3d
1352! Attention : la version precedente n'etait pas tres propre.
1353! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1354! pour obtenir le meme resultat.
1355         dtime=pdtphys
1356         radpas = NINT( 86400./dtime/nbapp_rad)
1357!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1358
1359         CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1360cIM begin
1361          print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
1362     $,ratqs(1,1)
1363cIM end
1364
1365
1366
1367!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1368c
1369C on remet le calendrier a zero
1370c
1371         IF (raz_date .eq. 1) THEN
1372           itau_phy = 0
1373         ENDIF
1374
1375cIM cf. AM 081204 BEG
1376         PRINT*,'cycle_diurne3 =',cycle_diurne
1377cIM cf. AM 081204 END
1378c
1379         CALL printflag( tabcntr0,radpas,ok_journe,
1380     ,                    ok_instan, ok_region )
1381c
1382         IF (ABS(dtime-pdtphys).GT.0.001) THEN
1383            WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1384     .                        pdtphys
1385            abort_message='Pas physique n est pas correct '
1386!           call abort_gcm(modname,abort_message,1)
1387            dtime=pdtphys
1388         ENDIF
1389         IF (nlon .NE. klon) THEN
1390            WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1391     .                      klon
1392            abort_message='nlon et klon ne sont pas coherents'
1393            call abort_gcm(modname,abort_message,1)
1394         ENDIF
1395         IF (nlev .NE. klev) THEN
1396            WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1397     .                       klev
1398            abort_message='nlev et klev ne sont pas coherents'
1399            call abort_gcm(modname,abort_message,1)
1400         ENDIF
1401c
1402         IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
1403           WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1404           WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1405           abort_message='Nbre d appels au rayonnement insuffisant'
1406           call abort_gcm(modname,abort_message,1)
1407         ENDIF
1408         WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1409         WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1410     .                   ok_cvl
1411c
1412cKE43
1413c Initialisation pour la convection de K.E. (sb):
1414         IF (iflag_con.GE.3) THEN
1415
1416         WRITE(lunout,*)"*** Convection de Kerry Emanuel 4.3  "
1417         WRITE(lunout,*)
1418     .      "On va utiliser le melange convectif des traceurs qui"
1419         WRITE(lunout,*)"est calcule dans convect4.3"
1420         WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1421
1422          DO i = 1, klon
1423           ema_cbmf(i) = 0.
1424           ema_pcb(i)  = 0.
1425           ema_pct(i)  = 0.
1426           ema_workcbmf(i) = 0.
1427          ENDDO
1428cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1429          DO i = 1, klon
1430           ibas_con(i) = 1
1431           itop_con(i) = 1
1432          ENDDO
1433cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1434c===============================================================================
1435cCR:04.12.07: initialisations poches froides
1436c Controle de ALE et ALP pour la fermeture convective (jyg)
1437          if (iflag_wake.eq.1) then
1438            CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
1439     s                  ,alp_bl_prescr, ale_bl_prescr)
1440c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1441c        print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1442          endif
1443
1444        do i = 1,klon
1445         Ale_bl(i)=0.
1446         Alp_bl(i)=0.
1447        enddo
1448
1449c================================================================================
1450
1451         ENDIF
1452
1453           DO i=1,klon
1454             rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1455           ENDDO
1456
1457c34EK
1458         IF (ok_orodr) THEN
1459
1460!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1461! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1462! justement quand ok_orodr = false.
1463! ce rugoro est utilise par la couche limite et fait double emploi
1464! avec les paramétrisations spécifiques de Francois Lott.
1465!           DO i=1,klon
1466!             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1467!           ENDDO
1468!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1469           IF (ok_strato) THEN
1470             CALL SUGWD_strato(klon,klev,paprs,pplay)
1471           ELSE
1472             CALL SUGWD(klon,klev,paprs,pplay)
1473           ENDIF
1474           
1475           DO i=1,klon
1476             zuthe(i)=0.
1477             zvthe(i)=0.
1478             if(zstd(i).gt.10.)then
1479               zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1480               zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1481             endif
1482           ENDDO
1483         ENDIF
1484c
1485c
1486         lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
1487         WRITE(lunout,*)'La frequence de lecture surface est de ',
1488     .                   lmt_pas
1489c
1490cIM 030306 END
1491
1492      capemaxcels = 't_max(X)'
1493      t2mincels = 't_min(X)'
1494      t2maxcels = 't_max(X)'
1495      tinst = 'inst(X)'
1496      tave = 'ave(X)'
1497cIM cf. AM 081204 BEG
1498      write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1499cIM cf. AM 081204 END
1500c
1501c=============================================================
1502c   Initialisation des sorties
1503c=============================================================
1504
1505#ifdef CPP_IOIPSL
1506
1507c$OMP MASTER
1508       call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta,
1509     &                        ctetaSTD,dtime,ok_veget,
1510     &                        type_ocean,iflag_pbl,ok_mensuel,ok_journe,
1511     &                        ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,
1512     &                        read_climoz, new_aod, aerosol_couple
1513     &                        )
1514c$OMP END MASTER
1515c$OMP BARRIER
1516
1517#ifdef histISCCP
1518#include "ini_histISCCP.h"
1519#endif
1520
1521#ifdef histNMC
1522#include "ini_histhfNMC.h"
1523#include "ini_histdayNMC.h"
1524#include "ini_histmthNMC.h"
1525#endif
1526
1527#include "ini_histday_seri.h"
1528
1529#include "ini_paramLMDZ_phy.h"
1530
1531#endif
1532
1533cIM 250308bad guide        ecrit_hf2mth = 30*1/ecrit_hf
1534         ecrit_hf2mth = ecrit_mth/ecrit_hf
1535
1536         ecrit_hf = ecrit_hf * un_jour
1537cIM
1538         IF(ecrit_day.LE.1.) THEN
1539          ecrit_day = ecrit_day * un_jour !en secondes
1540         ENDIF
1541cIM
1542         ecrit_mth = ecrit_mth * un_jour
1543         ecrit_ins = ecrit_ins * un_jour
1544         ecrit_reg = ecrit_reg * un_jour
1545         ecrit_tra = ecrit_tra * un_jour
1546         ecrit_ISCCP = ecrit_ISCCP * un_jour
1547         ecrit_LES = ecrit_LES * un_jour
1548c
1549         PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1550     .   ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1551     .   ecrit_hf2mth
1552cIM 030306 END
1553
1554
1555cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1556      date0 = jD_ref
1557      WRITE(*,*) 'physiq date0 : ',date0
1558c
1559c
1560c
1561c Prescrire l'ozone dans l'atmosphere
1562c
1563c
1564cc         DO i = 1, klon
1565cc         DO k = 1, klev
1566cc            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1567cc         ENDDO
1568cc         ENDDO
1569c
1570      IF (config_inca /= 'none') THEN
1571#ifdef INCA
1572         CALL VTe(VTphysiq)
1573         CALL VTb(VTinca)
1574!         iii = MOD(NINT(xjour),360)
1575!         calday = REAL(iii) + jH_cur
1576         calday = REAL(days_elapsed) + jH_cur
1577         WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1578
1579         CALL chemini(
1580     $                   rg,
1581     $                   ra,
1582     $                   airephy,
1583     $                   rlat,
1584     $                   rlon,
1585     $                   presnivs,
1586     $                   calday,
1587     $                   klon,
1588     $                   nqtot,
1589     $                   pdtphys,
1590     $                   annee_ref,
1591     $                   day_ref,
1592     $                   itau_phy)
1593
1594         CALL VTe(VTinca)
1595         CALL VTb(VTphysiq)
1596#endif
1597      END IF
1598c
1599c
1600!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1601! Nouvelle initialisation pour le rayonnement RRTM
1602!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1603
1604      call iniradia(klon,klev,paprs(1,1:klev+1))
1605
1606C$omp single
1607      if (read_climoz >= 1) then
1608         call open_climoz(ncid_climoz, press_climoz)
1609      END IF
1610C$omp end single
1611      ENDIF
1612!
1613!   ****************     Fin  de   IF ( debut  )   ***************
1614!
1615!
1616! Incrementer le compteur de la physique
1617!
1618      itap   = itap + 1
1619!
1620! Update fraction of the sub-surfaces (pctsrf) and
1621! initialize, where a new fraction has appeared, all variables depending
1622! on the surface fraction.
1623!
1624      CALL change_srf_frac(itap, dtime, days_elapsed+1,
1625     *     pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
1626
1627! Tendances bidons pour les processus qui n'affectent pas certaines
1628! variables.
1629      du0(:,:)=0.
1630      dv0(:,:)=0.
1631      dq0(:,:)=0.
1632      dql0(:,:)=0.
1633c
1634c Mettre a zero des variables de sortie (pour securite)
1635c
1636      DO i = 1, klon
1637         d_ps(i) = 0.0
1638      ENDDO
1639      DO k = 1, klev
1640      DO i = 1, klon
1641         d_t(i,k) = 0.0
1642         d_u(i,k) = 0.0
1643         d_v(i,k) = 0.0
1644      ENDDO
1645      ENDDO
1646      DO iq = 1, nqtot
1647      DO k = 1, klev
1648      DO i = 1, klon
1649         d_qx(i,k,iq) = 0.0
1650      ENDDO
1651      ENDDO
1652      ENDDO
1653      da(:,:)=0.
1654      mp(:,:)=0.
1655      phi(:,:,:)=0.
1656c
1657c Ne pas affecter les valeurs entrees de u, v, h, et q
1658c
1659      DO k = 1, klev
1660      DO i = 1, klon
1661         t_seri(i,k)  = t(i,k)
1662         u_seri(i,k)  = u(i,k)
1663         v_seri(i,k)  = v(i,k)
1664         q_seri(i,k)  = qx(i,k,ivap)
1665         ql_seri(i,k) = qx(i,k,iliq)
1666         qs_seri(i,k) = 0.
1667      ENDDO
1668      ENDDO
1669      IF (nqtot.GE.3) THEN
1670      DO iq = 3, nqtot
1671      DO  k = 1, klev
1672      DO  i = 1, klon
1673         tr_seri(i,k,iq-2) = qx(i,k,iq)
1674      ENDDO
1675      ENDDO
1676      ENDDO
1677      ELSE
1678      DO k = 1, klev
1679      DO i = 1, klon
1680         tr_seri(i,k,1) = 0.0
1681      ENDDO
1682      ENDDO
1683      ENDIF
1684C
1685      DO i = 1, klon
1686        ztsol(i) = 0.
1687      ENDDO
1688      DO nsrf = 1, nbsrf
1689        DO i = 1, klon
1690          ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1691        ENDDO
1692      ENDDO
1693cIM
1694      IF (ip_ebil_phy.ge.1) THEN
1695        ztit='after dynamic'
1696        CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1697     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1698     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1699C     Comme les tendances de la physique sont ajoute dans la dynamique,
1700C     on devrait avoir que la variation d'entalpie par la dynamique
1701C     est egale a la variation de la physique au pas de temps precedent.
1702C     Donc la somme de ces 2 variations devrait etre nulle.
1703        call diagphy(airephy,ztit,ip_ebil_phy
1704     e      , zero_v, zero_v, zero_v, zero_v, zero_v
1705     e      , zero_v, zero_v, zero_v, ztsol
1706     e      , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1707     s      , fs_bound, fq_bound )
1708      END IF
1709
1710c Diagnostiquer la tendance dynamique
1711c
1712      IF (ancien_ok) THEN
1713         DO k = 1, klev
1714         DO i = 1, klon
1715            d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1716            d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1717            d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1718            d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1719         ENDDO
1720         ENDDO
1721      ELSE
1722         DO k = 1, klev
1723         DO i = 1, klon
1724            d_u_dyn(i,k) = 0.0
1725            d_v_dyn(i,k) = 0.0
1726            d_t_dyn(i,k) = 0.0
1727            d_q_dyn(i,k) = 0.0
1728         ENDDO
1729         ENDDO
1730         ancien_ok = .TRUE.
1731      ENDIF
1732c
1733c Ajouter le geopotentiel du sol:
1734c
1735      DO k = 1, klev
1736      DO i = 1, klon
1737         zphi(i,k) = pphi(i,k) + pphis(i)
1738      ENDDO
1739      ENDDO
1740c
1741c Verifier les temperatures
1742c
1743cIM BEG
1744      IF (check) THEN
1745       amn=MIN(ftsol(1,is_ter),1000.)
1746       amx=MAX(ftsol(1,is_ter),-1000.)
1747       DO i=2, klon
1748        amn=MIN(ftsol(i,is_ter),amn)
1749        amx=MAX(ftsol(i,is_ter),amx)
1750       ENDDO
1751c
1752       PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1753      ENDIF !(check) THEN
1754cIM END
1755c
1756      CALL hgardfou(t_seri,ftsol,'debutphy')
1757c
1758cIM BEG
1759      IF (check) THEN
1760       amn=MIN(ftsol(1,is_ter),1000.)
1761       amx=MAX(ftsol(1,is_ter),-1000.)
1762       DO i=2, klon
1763        amn=MIN(ftsol(i,is_ter),amn)
1764        amx=MAX(ftsol(i,is_ter),amx)
1765       ENDDO
1766c
1767       PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1768      ENDIF !(check) THEN
1769cIM END
1770c
1771c Mettre en action les conditions aux limites (albedo, sst, etc.).
1772c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1773c
1774      if (read_climoz >= 1) then
1775C        Ozone from a file
1776!        Update required ozone index:
1777         ro3i = int((days_elapsed + jh_cur - jh_1jan)
1778     $        / ioget_year_len(year_cur) * 360.) + 1
1779         if (ro3i == 361) ro3i = 360
1780C        (This should never occur, except perhaps because of roundup
1781C        error. See documentation.)
1782         if (ro3i /= co3i) then
1783C           Update ozone field:
1784            if (read_climoz == 1) then
1785               call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
1786     $              press_in_edg=press_climoz, paprs=paprs, v3=wo)
1787            else
1788C              read_climoz == 2
1789               call regr_pr_av(ncid_climoz,
1790     $              (/"tro3         ", "tro3_daylight"/),
1791     $              julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
1792     $              v3=wo)
1793            end if
1794!           Convert from mole fraction of ozone to column density of ozone in a
1795!           cell, in kDU:
1796            forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
1797     $           * rmo3 / rmd * zmasse / dobson_u / 1e3
1798C           (By regridding ozone values for LMDZ only once every 360th of
1799C           year, we have already neglected the variation of pressure in one
1800C           360th of year. So do not recompute "wo" at each time step even if
1801C           "zmasse" changes a little.)
1802            co3i = ro3i
1803         end if
1804      elseif (MOD(itap-1,lmt_pas) == 0) THEN
1805C        Once per day, update ozone from Royer:
1806         wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
1807      ENDIF
1808c
1809c Re-evaporer l'eau liquide nuageuse
1810c
1811      DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
1812      DO i = 1, klon
1813         zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1814c        zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1815         zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1816         zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1817         zb = MAX(0.0,ql_seri(i,k))
1818         za = - MAX(0.0,ql_seri(i,k))
1819     .                  * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
1820         t_seri(i,k) = t_seri(i,k) + za
1821         q_seri(i,k) = q_seri(i,k) + zb
1822         ql_seri(i,k) = 0.0
1823         d_t_eva(i,k) = za
1824         d_q_eva(i,k) = zb
1825      ENDDO
1826      ENDDO
1827cIM
1828      IF (ip_ebil_phy.ge.2) THEN
1829        ztit='after reevap'
1830        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1831     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1832     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1833         call diagphy(airephy,ztit,ip_ebil_phy
1834     e      , zero_v, zero_v, zero_v, zero_v, zero_v
1835     e      , zero_v, zero_v, zero_v, ztsol
1836     e      , d_h_vcol, d_qt, d_ec
1837     s      , fs_bound, fq_bound )
1838C
1839      END IF
1840
1841c
1842c=========================================================================
1843! Calculs de l'orbite.
1844! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1845! doit donc etre placé avant radlwsw et pbl_surface
1846
1847! calcul selon la routine utilisee pour les planetes
1848      if (new_orbit) then
1849        call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1850        day_since_equinox = (jD_cur + jH_cur) - jD_eq
1851!        day_since_equinox = (jD_cur) - jD_eq
1852        call solarlong(day_since_equinox, zlongi, dist)
1853      else     
1854! calcul selon la routine utilisee pour l'AR4
1855!   choix entre calcul de la longitude solaire vraie ou valeur fixee a
1856!   solarlong0
1857        if (solarlong0<-999.) then
1858           CALL orbite(REAL(days_elapsed+1),zlongi,dist)
1859        else
1860           zlongi=solarlong0  ! longitude solaire vraie
1861           dist=1.            ! distance au soleil / moyenne
1862        endif
1863      endif
1864      if(prt_level.ge.1)                                                &
1865     &    write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1866
1867!  Avec ou sans cycle diurne
1868      IF (cycle_diurne) THEN
1869        zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1870        CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
1871      ELSE
1872        CALL angle(zlongi, rlat, fract, rmu0)
1873      ENDIF
1874
1875      if (mydebug) then
1876        call writefield_phy('u_seri',u_seri,llm)
1877        call writefield_phy('v_seri',v_seri,llm)
1878        call writefield_phy('t_seri',t_seri,llm)
1879        call writefield_phy('q_seri',q_seri,llm)
1880      endif
1881
1882ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1883c Appel au pbl_surface : Planetary Boudary Layer et Surface
1884c Cela implique tous les interactions des sous-surfaces et la partie diffusion
1885c turbulent du couche limit.
1886c
1887c Certains varibales de sorties de pbl_surface sont utiliser que pour
1888c ecriture des fihiers hist_XXXX.nc, ces sont :
1889c   qsol,      zq2m,      s_pblh,  s_lcl,
1890c   s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
1891c   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
1892c   zxrugs,    zu10m,     zv10m,   fder,
1893c   zxqsurf,   rh2m,      zxfluxu, zxfluxv,
1894c   frugs,     agesno,    fsollw,  fsolsw,
1895c   d_ts,      fevap,     fluxlat, t2m,
1896c   wfbils,    wfbilo,    fluxt,   fluxu, fluxv,
1897c
1898c Certains ne sont pas utiliser du tout :
1899c   dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1900c
1901
1902      CALL pbl_surface(
1903     e     dtime,     date0,     itap,    days_elapsed+1,
1904     e     debut,     lafin,
1905     e     rlon,      rlat,      rugoro,  rmu0,     
1906     e     rain_fall, snow_fall, solsw,   sollw,   
1907     e     t_seri,    q_seri,    u_seri,  v_seri,   
1908     e     pplay,     paprs,     pctsrf,           
1909     +     ftsol,     falb1,     falb2,   u10m,   v10m,
1910     s     sollwdown, cdragh,    cdragm,  u1,    v1,
1911     s     albsol1,   albsol2,   sens,    evap, 
1912     s     zxtsol,    zxfluxlat, zt2m,    qsat2m,
1913     s     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
1914     s     coefh,     slab_wfbils,               
1915     d     qsol,      zq2m,      s_pblh,  s_lcl,
1916     d     s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
1917     d     s_therm,   s_trmb1,   s_trmb2, s_trmb3,
1918     d     zxrugs,    zu10m,     zv10m,   fder,
1919     d     zxqsurf,   rh2m,      zxfluxu, zxfluxv,
1920     d     frugs,     agesno,    fsollw,  fsolsw,
1921     d     d_ts,      fevap,     fluxlat, t2m,
1922     d     wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
1923     -     dsens,     devap,     zxsnow,
1924     -     zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
1925
1926
1927!-----------------------------------------------------------------------------------------
1928! ajout des tendances de la diffusion turbulente
1929      CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
1930!-----------------------------------------------------------------------------------------
1931
1932      if (mydebug) then
1933        call writefield_phy('u_seri',u_seri,llm)
1934        call writefield_phy('v_seri',v_seri,llm)
1935        call writefield_phy('t_seri',t_seri,llm)
1936        call writefield_phy('q_seri',q_seri,llm)
1937      endif
1938
1939
1940      IF (ip_ebil_phy.ge.2) THEN
1941        ztit='after surface_main'
1942        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1943     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1944     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1945         call diagphy(airephy,ztit,ip_ebil_phy
1946     e      , zero_v, zero_v, zero_v, zero_v, sens
1947     e      , evap  , zero_v, zero_v, ztsol
1948     e      , d_h_vcol, d_qt, d_ec
1949     s      , fs_bound, fq_bound )
1950      END IF
1951
1952c =================================================================== c
1953c   Calcul de Qsat
1954
1955      DO k = 1, klev
1956      DO i = 1, klon
1957         zx_t = t_seri(i,k)
1958         IF (thermcep) THEN
1959            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1960            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1961            zx_qs  = MIN(0.5,zx_qs)
1962            zcor   = 1./(1.-retv*zx_qs)
1963            zx_qs  = zx_qs*zcor
1964         ELSE
1965           IF (zx_t.LT.t_coup) THEN
1966              zx_qs = qsats(zx_t)/pplay(i,k)
1967           ELSE
1968              zx_qs = qsatl(zx_t)/pplay(i,k)
1969           ENDIF
1970         ENDIF
1971         zqsat(i,k)=zx_qs
1972      ENDDO
1973      ENDDO
1974
1975      if (prt_level.ge.1) then
1976      write(lunout,*) 'L   qsat (g/kg) avant clouds_gno'
1977      write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1978      endif
1979c
1980c Appeler la convection (au choix)
1981c
1982      DO k = 1, klev
1983      DO i = 1, klon
1984         conv_q(i,k) = d_q_dyn(i,k)
1985     .               + d_q_vdf(i,k)/dtime
1986         conv_t(i,k) = d_t_dyn(i,k)
1987     .               + d_t_vdf(i,k)/dtime
1988      ENDDO
1989      ENDDO
1990      IF (check) THEN
1991         za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1992         WRITE(lunout,*) "avantcon=", za
1993      ENDIF
1994      zx_ajustq = .FALSE.
1995      IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1996      IF (zx_ajustq) THEN
1997         DO i = 1, klon
1998            z_avant(i) = 0.0
1999         ENDDO
2000         DO k = 1, klev
2001         DO i = 1, klon
2002            z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
2003     .                        *(paprs(i,k)-paprs(i,k+1))/RG
2004         ENDDO
2005         ENDDO
2006      ENDIF
2007
2008c Calcule de vitesse verticale a partir de flux de masse verticale
2009      DO k = 1, klev
2010         DO i = 1, klon
2011            omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
2012         END DO
2013      END DO
2014      if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ',
2015     $     omega(igout, :)
2016
2017      IF (iflag_con.EQ.1) THEN
2018        abort_message ='reactiver le call conlmd dans physiq.F'
2019        CALL abort_gcm (modname,abort_message,1)
2020c     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2021c    .             d_t_con, d_q_con,
2022c    .             rain_con, snow_con, ibas_con, itop_con)
2023      ELSE IF (iflag_con.EQ.2) THEN
2024      CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
2025     e            conv_t, conv_q, -evap, omega,
2026     s            d_t_con, d_q_con, rain_con, snow_con,
2027     s            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2028     s            kcbot, kctop, kdtop, pmflxr, pmflxs)
2029      d_u_con = 0.
2030      d_v_con = 0.
2031
2032      WHERE (rain_con < 0.) rain_con = 0.
2033      WHERE (snow_con < 0.) snow_con = 0.
2034      DO i = 1, klon
2035         ibas_con(i) = klev+1 - kcbot(i)
2036         itop_con(i) = klev+1 - kctop(i)
2037      ENDDO
2038      ELSE IF (iflag_con.GE.3) THEN
2039c nb of tracers for the KE convection:
2040c MAF la partie traceurs est faite dans phytrac
2041c on met ntra=1 pour limiter les appels mais on peut
2042c supprimer les calculs / ftra.
2043              ntra = 1
2044
2045c=====================================================================================
2046cajout pour la parametrisation des poches froides:
2047ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2048      do k=1,klev
2049            do i=1,klon
2050             if (iflag_wake.eq.1) then
2051             t_wake(i,k) = t_seri(i,k)
2052     .           +(1-wake_s(i))*wake_deltat(i,k)
2053             q_wake(i,k) = q_seri(i,k)
2054     .           +(1-wake_s(i))*wake_deltaq(i,k)
2055             t_undi(i,k) = t_seri(i,k)
2056     .           -wake_s(i)*wake_deltat(i,k)
2057             q_undi(i,k) = q_seri(i,k)
2058     .           -wake_s(i)*wake_deltaq(i,k)
2059             else
2060             t_wake(i,k) = t_seri(i,k)
2061             q_wake(i,k) = q_seri(i,k)
2062             t_undi(i,k) = t_seri(i,k)
2063             q_undi(i,k) = q_seri(i,k)
2064             endif
2065            enddo
2066         enddo
2067     
2068cc--   Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2069cc--    pour le soulevement des particules dans le modele convectif
2070c
2071      do i = 1,klon
2072        ALE(i) = 0.
2073        ALP(i) = 0.
2074      enddo
2075c
2076ccalcul de ale_wake et alp_wake
2077       if (iflag_wake.eq.1) then
2078         if (itap .le. it_wape_prescr) then
2079          do i = 1,klon
2080           ale_wake(i) = wape_prescr
2081           alp_wake(i) = fip_prescr
2082          enddo
2083         else
2084          do i = 1,klon
2085cjyg  ALE=WAPE au lieu de ALE = 1/2 Cstar**2
2086ccc           ale_wake(i) = 0.5*wake_cstar(i)**2
2087           ale_wake(i) = wake_pe(i)
2088           alp_wake(i) = wake_fip(i)
2089          enddo
2090         endif
2091       else
2092         do i = 1,klon
2093           ale_wake(i) = 0.
2094           alp_wake(i) = 0.
2095         enddo
2096       endif
2097ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2098cdans le thermique sinon
2099       if (iflag_coupl.eq.0) then
2100          if (debut.and.prt_level.gt.9)
2101     $                     WRITE(lunout,*)'ALE et ALP imposes'
2102          do i = 1,klon
2103con ne couple que ale
2104c           ALE(i) = max(ale_wake(i),Ale_bl(i))
2105            ALE(i) = max(ale_wake(i),ale_bl_prescr)
2106con ne couple que alp
2107c           ALP(i) = alp_wake(i) + Alp_bl(i)
2108            ALP(i) = alp_wake(i) + alp_bl_prescr
2109          enddo
2110       else
2111         IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2112!         do i = 1,klon
2113!             ALE(i) = max(ale_wake(i),Ale_bl(i))
2114! avant        ALP(i) = alp_wake(i) + Alp_bl(i)
2115!             ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2116!         write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2117!         write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2118!         enddo
2119
2120!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2121! Modif FH 2010/04/27. Sans doute temporaire.
2122! Deux options pour le alp_offset : constant si >Â 0 ou proportionnel Ãa
2123! w si <0
2124!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2125       do i = 1,klon
2126          ALE(i) = max(ale_wake(i),Ale_bl(i))
2127          if (alp_offset>=0.) then
2128            ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2129          else
2130            ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2131            if (alp(i)<0.) then
2132               print*,'ALP ',alp(i),alp_wake(i)
2133     s         ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2134            endif
2135          endif
2136       enddo
2137!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2138
2139       endif
2140       do i=1,klon
2141          if (alp(i)>alp_max) then
2142             IF(prt_level>9)WRITE(lunout,*)                             &
2143     &       'WARNING SUPER ALP (seuil=',alp_max,
2144     ,       '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2145             alp(i)=alp_max
2146          endif
2147          if (ale(i)>ale_max) then
2148             IF(prt_level>9)WRITE(lunout,*)                             &
2149     &       'WARNING SUPER ALE (seuil=',ale_max,
2150     ,       '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2151             ale(i)=ale_max
2152          endif
2153       enddo
2154
2155cfin calcul ale et alp
2156c=================================================================================================
2157
2158
2159c sb, oct02:
2160c Schema de convection modularise et vectorise:
2161c (driver commun aux versions 3 et 4)
2162c
2163          IF (ok_cvl) THEN ! new driver for convectL
2164
2165          CALL concvl (iflag_con,iflag_clos,
2166     .        dtime,paprs,pplay,t_undi,q_undi,
2167     .        t_wake,q_wake,wake_s,
2168     .        u_seri,v_seri,tr_seri,nbtr,
2169     .        ALE,ALP,
2170     .        ema_work1,ema_work2,
2171     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2172     .        rain_con, snow_con, ibas_con, itop_con, sigd,
2173     .        upwd,dnwd,dnwd0,
2174     .        Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2175     .        pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2176     .        pmflxr,pmflxs,da,phi,mp,
2177     .        ftd,fqd,lalim_conv,wght_th)
2178
2179cIM begin
2180c       print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2181c    .dnwd0(1,1),ftd(1,1),fqd(1,1)
2182cIM end
2183cIM cf. FH
2184              clwcon0=qcondc
2185              pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2186
2187              do i = 1, klon
2188                if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2189              enddo
2190
2191          ELSE ! ok_cvl
2192c MAF conema3 ne contient pas les traceurs
2193          CALL conema3 (dtime,
2194     .        paprs,pplay,t_seri,q_seri,
2195     .        u_seri,v_seri,tr_seri,ntra,
2196     .        ema_work1,ema_work2,
2197     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2198     .        rain_con, snow_con, ibas_con, itop_con,
2199     .        upwd,dnwd,dnwd0,bas,top,
2200     .        Ma,cape,tvp,rflag,
2201     .        pbase
2202     .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2203     .        ,clwcon0)
2204
2205          ENDIF ! ok_cvl
2206
2207c
2208c Correction precip
2209          rain_con = rain_con * cvl_corr
2210          snow_con = snow_con * cvl_corr
2211c
2212
2213           IF (.NOT. ok_gust) THEN
2214           do i = 1, klon
2215            wd(i)=0.0
2216           enddo
2217           ENDIF
2218
2219c =================================================================== c
2220c Calcul des proprietes des nuages convectifs
2221c
2222
2223c   calcul des proprietes des nuages convectifs
2224             clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2225             call clouds_gno
2226     s       (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2227
2228c =================================================================== c
2229
2230          DO i = 1, klon
2231            ema_pcb(i)  = pbase(i)
2232          ENDDO
2233          DO i = 1, klon
2234
2235! L'idicage de itop_con peut cacher un pb potentiel
2236! FH sous la dictee de JYG, CR
2237            ema_pct(i)  = paprs(i,itop_con(i)+1)
2238
2239            if (itop_con(i).gt.klev-3) then
2240              if(prt_level >= 9) then
2241                write(lunout,*)'La convection monte trop haut '
2242                write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2243              endif
2244            endif
2245          ENDDO
2246          DO i = 1, klon
2247            ema_cbmf(i) = ema_workcbmf(i)
2248          ENDDO     
2249      ELSE IF (iflag_con.eq.0) THEN
2250          write(lunout,*) 'On n appelle pas la convection'
2251          clwcon0=0.
2252          rnebcon0=0.
2253          d_t_con=0.
2254          d_q_con=0.
2255          d_u_con=0.
2256          d_v_con=0.
2257          rain_con=0.
2258          snow_con=0.
2259          bas=1
2260          top=1
2261      ELSE
2262          WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2263          CALL abort
2264      ENDIF
2265
2266c     CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2267c    .              d_u_con, d_v_con)
2268
2269!-----------------------------------------------------------------------------------------
2270! ajout des tendances de la diffusion turbulente
2271      CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2272!-----------------------------------------------------------------------------------------
2273
2274      if (mydebug) then
2275        call writefield_phy('u_seri',u_seri,llm)
2276        call writefield_phy('v_seri',v_seri,llm)
2277        call writefield_phy('t_seri',t_seri,llm)
2278        call writefield_phy('q_seri',q_seri,llm)
2279      endif
2280
2281cIM
2282      IF (ip_ebil_phy.ge.2) THEN
2283        ztit='after convect'
2284        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2285     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2286     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2287         call diagphy(airephy,ztit,ip_ebil_phy
2288     e      , zero_v, zero_v, zero_v, zero_v, zero_v
2289     e      , zero_v, rain_con, snow_con, ztsol
2290     e      , d_h_vcol, d_qt, d_ec
2291     s      , fs_bound, fq_bound )
2292      END IF
2293C
2294      IF (check) THEN
2295          za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2296          WRITE(lunout,*)"aprescon=", za
2297          zx_t = 0.0
2298          za = 0.0
2299          DO i = 1, klon
2300            za = za + airephy(i)/REAL(klon)
2301            zx_t = zx_t + (rain_con(i)+
2302     .                   snow_con(i))*airephy(i)/REAL(klon)
2303          ENDDO
2304          zx_t = zx_t/za*dtime
2305          WRITE(lunout,*)"Precip=", zx_t
2306      ENDIF
2307      IF (zx_ajustq) THEN
2308          DO i = 1, klon
2309            z_apres(i) = 0.0
2310          ENDDO
2311          DO k = 1, klev
2312            DO i = 1, klon
2313              z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2314     .            *(paprs(i,k)-paprs(i,k+1))/RG
2315            ENDDO
2316          ENDDO
2317          DO i = 1, klon
2318            z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2319     .          /z_apres(i)
2320          ENDDO
2321          DO k = 1, klev
2322            DO i = 1, klon
2323              IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2324     .            z_factor(i).LT.(1.0-1.0E-08)) THEN
2325                  q_seri(i,k) = q_seri(i,k) * z_factor(i)
2326              ENDIF
2327            ENDDO
2328          ENDDO
2329      ENDIF
2330      zx_ajustq=.FALSE.
2331
2332c
2333c=============================================================================
2334cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2335cpour la couche limite diffuse pour l instant
2336c
2337      if (iflag_wake.eq.1) then
2338      DO k=1,klev
2339        DO i=1,klon
2340          dt_dwn(i,k)  = ftd(i,k)
2341          wdt_PBL(i,k) = 0.
2342          dq_dwn(i,k)  = fqd(i,k)
2343          wdq_PBL(i,k) = 0.
2344          M_dwn(i,k)   = dnwd0(i,k)
2345          M_up(i,k)    = upwd(i,k)
2346          dt_a(i,k)    = d_t_con(i,k)/dtime - ftd(i,k)
2347          udt_PBL(i,k) = 0.
2348          dq_a(i,k)    = d_q_con(i,k)/dtime - fqd(i,k)
2349          udq_PBL(i,k) = 0.
2350        ENDDO
2351      ENDDO
2352c
2353ccalcul caracteristiques de la poche froide
2354      call calWAKE (paprs,pplay,dtime
2355     :               ,t_seri,q_seri,omega
2356     :               ,dt_dwn,dq_dwn,M_dwn,M_up
2357     :               ,dt_a,dq_a,sigd
2358     :               ,wdt_PBL,wdq_PBL
2359     :               ,udt_PBL,udq_PBL
2360     o               ,wake_deltat,wake_deltaq,wake_dth
2361     o               ,wake_h,wake_s,wake_dens
2362     o               ,wake_pe,wake_fip,wake_gfl
2363     o               ,dt_wake,dq_wake
2364     o               ,wake_k, t_undi,q_undi
2365     o               ,wake_omgbdth,wake_dp_omgb
2366     o               ,wake_dtKE,wake_dqKE
2367     o               ,wake_dtPBL,wake_dqPBL
2368     o               ,wake_omg,wake_dp_deltomg
2369     o               ,wake_spread,wake_Cstar,wake_d_deltat_gw
2370     o               ,wake_ddeltat,wake_ddeltaq)
2371c
2372!-----------------------------------------------------------------------------------------
2373! ajout des tendances des poches froides
2374! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2375! coherent avec les autres d_t_...
2376      d_t_wake(:,:)=dt_wake(:,:)*dtime
2377      d_q_wake(:,:)=dq_wake(:,:)*dtime
2378      CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2379!-----------------------------------------------------------------------------------------
2380
2381      endif
2382c
2383c===================================================================
2384cJYG
2385      IF (ip_ebil_phy.ge.2) THEN
2386        ztit='after wake'
2387        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2388     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2389     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2390        call diagphy(airephy,ztit,ip_ebil_phy
2391     e      , zero_v, zero_v, zero_v, zero_v, zero_v
2392     e      , zero_v, zero_v, zero_v, ztsol
2393     e      , d_h_vcol, d_qt, d_ec
2394     s      , fs_bound, fq_bound )
2395      END IF
2396
2397c      print*,'apres callwake iflag_cldcon=', iflag_cldcon
2398c
2399c===================================================================
2400c Convection seche (thermiques ou ajustement)
2401c===================================================================
2402c
2403       call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2404     s ,seuil_inversion,weak_inversion,dthmin)
2405
2406
2407
2408      d_t_ajsb(:,:)=0.
2409      d_q_ajsb(:,:)=0.
2410      d_t_ajs(:,:)=0.
2411      d_u_ajs(:,:)=0.
2412      d_v_ajs(:,:)=0.
2413      d_q_ajs(:,:)=0.
2414      clwcon0th(:,:)=0.
2415c
2416c      fm_therm(:,:)=0.
2417c      entr_therm(:,:)=0.
2418c      detr_therm(:,:)=0.
2419c
2420      IF(prt_level>9)WRITE(lunout,*)
2421     .    'AVANT LA CONVECTION SECHE , iflag_thermals='
2422     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
2423      if(iflag_thermals.lt.0) then
2424c  Rien
2425c  ====
2426         IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2427
2428
2429      else
2430
2431c  Thermiques
2432c  ==========
2433         IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2434     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
2435
2436
2437         if (iflag_thermals.gt.1) then
2438         call calltherm(pdtphys
2439     s      ,pplay,paprs,pphi,weak_inversion
2440     s      ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2441     s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2442     s      ,fm_therm,entr_therm,detr_therm
2443     s      ,zqasc,clwcon0th,lmax_th,ratqscth
2444     s      ,ratqsdiff,zqsatth
2445con rajoute ale et alp, et les caracteristiques de la couche alim
2446     s      ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca
2447     s      ,ztv,zpspsk,ztla,zthl)
2448
2449! ----------------------------------------------------------------------
2450! Transport de la TKE par les panaches thermiques.
2451! FH : 2010/02/01
2452      if (iflag_pbl.eq.10) then
2453      call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
2454     s           rg,paprs,pbl_tke)
2455      endif
2456! ----------------------------------------------------------------------
2457
2458         endif
2459
2460
2461
2462c  Ajustement sec
2463c  ==============
2464
2465! Dans le cas où on active les thermiques, on fait partir l'ajustement
2466! a partir du sommet des thermiques.
2467! Dans le cas contraire, on demarre au niveau 1.
2468
2469         if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2470
2471         if(iflag_thermals.eq.0) then
2472            IF(prt_level>9)WRITE(lunout,*)'ajsec'
2473            limbas(:)=1
2474         else
2475            limbas(:)=lmax_th(:)
2476         endif
2477 
2478! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2479! pour des test de convergence numerique.
2480! Le nouveau ajsec est a priori mieux, meme pour le cas
2481! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2482! non nulles numeriquement pour des mailles non concernees.
2483
2484         if (iflag_thermals.eq.0) then
2485            CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2486     s      , d_t_ajsb, d_q_ajsb)
2487         else
2488            CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2489     s      , d_t_ajsb, d_q_ajsb)
2490         endif
2491
2492!-----------------------------------------------------------------------------------------
2493! ajout des tendances de l'ajustement sec ou des thermiques
2494      CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2495         d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2496         d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2497
2498!-----------------------------------------------------------------------------------------
2499
2500         endif
2501
2502      endif
2503c
2504c===================================================================
2505cIM
2506      IF (ip_ebil_phy.ge.2) THEN
2507        ztit='after dry_adjust'
2508        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2509     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2510     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2511        call diagphy(airephy,ztit,ip_ebil_phy
2512     e      , zero_v, zero_v, zero_v, zero_v, zero_v
2513     e      , zero_v, zero_v, zero_v, ztsol
2514     e      , d_h_vcol, d_qt, d_ec
2515     s      , fs_bound, fq_bound )
2516      END IF
2517
2518
2519c-------------------------------------------------------------------------
2520c  Caclul des ratqs
2521c-------------------------------------------------------------------------
2522
2523c      print*,'calcul des ratqs'
2524c   ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2525c   ----------------
2526c   on ecrase le tableau ratqsc calcule par clouds_gno
2527      if (iflag_cldcon.eq.1) then
2528         do k=1,klev
2529         do i=1,klon
2530            if(ptconv(i,k)) then
2531              ratqsc(i,k)=ratqsbas
2532     s        +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2533            else
2534               ratqsc(i,k)=0.
2535            endif
2536         enddo
2537         enddo
2538
2539c-----------------------------------------------------------------------
2540c  par nversion de la fonction log normale
2541c-----------------------------------------------------------------------
2542      else if (iflag_cldcon.eq.4) then
2543         ptconvth(:,:)=.false.
2544         ratqsc(:,:)=0.
2545         if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2546         call clouds_gno
2547     s   (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2548         if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2549
2550c-----------------------------------------------------------------------
2551c   par calcul direct de l'ecart-type
2552c-----------------------------------------------------------------------
2553
2554      else if (iflag_cldcon>=5) then
2555         wmax_th(:)=0.
2556         zmax_th(:)=0.
2557         do k=1,klev
2558            do i=1,klon
2559               wmax_th(i)=max(wmax_th(i),zw2(i,k))
2560               if (detr_therm(i,k).gt.0.) zmax_th(i)=pphi(i,k)/rg
2561            enddo
2562         enddo
2563         tau_overturning_th(:)=zmax_th(:)/max(0.5*wmax_th(:),0.1)
2564         if(prt_level.ge.9)
2565     &     write(lunout,*)'TAU TH OK ',
2566     &     tau_overturning_th(1),detr_therm(1,3)
2567
2568c On impose que l'air autour de la fraction couverte par le thermique
2569c plus son air detraine durant tau_overturning_th soit superieur
2570c a 0.1 q_seri
2571         zz=0.1
2572         do k=1,klev
2573            do i=1,klon
2574               lambda_th(i,k)=0.5*(fraca(i,k)+fraca(i,k+1))+
2575     s         tau_overturning_th(i)*detr_therm(i,k)
2576     s         *rg/(paprs(i,k)-paprs(i,k+1))
2577               znum=(1.-zz)*q_seri(i,k)
2578               zden=zqasc(i,k)-zz*q_seri(i,k)
2579               if (znum-lambda_th(i,k)*zden<0.) lambda_th(i,k)=znum/zden
2580               lambda_th(i,k)=min(lambda_th(i,k),0.9)
2581            enddo
2582         enddo
2583
2584         if(iflag_cldcon==5) then
2585            do k=1,klev
2586               do i=1,klon
2587                  ratqsc(i,k)=sqrt(lambda_th(i,k)/(1.-lambda_th(i,k)))*
2588     s            abs((zqasc(i,k)-q_seri(i,k))/q_seri(i,k))
2589               enddo
2590            enddo
2591         else if(iflag_cldcon==6) then
2592            do k=1,klev
2593               do i=1,klon
2594                  ratqsc(i,k)=sqrt(lambda_th(i,k))*
2595     s            (zqasc(i,k)-q_seri(i,k))/q_seri(i,k)
2596               enddo
2597            enddo
2598         endif
2599
2600      endif
2601
2602c   ratqs stables
2603c   -------------
2604
2605      if (iflag_ratqs.eq.0) then
2606
2607! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2608         do k=1,klev
2609            do i=1, klon
2610               ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2611     s         min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2612            enddo
2613         enddo
2614
2615! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2616! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2617! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2618! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2619! Il s'agit de differents tests dans la phase de reglage du modele
2620! avec thermiques.
2621
2622      else if (iflag_ratqs.eq.1) then
2623
2624         do k=1,klev
2625            do i=1, klon
2626               if (pplay(i,k).ge.60000.) then
2627                  ratqss(i,k)=ratqsbas
2628               else if ((pplay(i,k).ge.30000.).and.
2629     s            (pplay(i,k).lt.60000.)) then
2630                  ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2631     s            (60000.-pplay(i,k))/(60000.-30000.)
2632               else
2633                  ratqss(i,k)=ratqshaut
2634               endif
2635            enddo
2636         enddo
2637
2638      else
2639
2640         do k=1,klev
2641            do i=1, klon
2642               if (pplay(i,k).ge.60000.) then
2643                  ratqss(i,k)=ratqsbas
2644     s            *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2645               else if ((pplay(i,k).ge.30000.).and.
2646     s             (pplay(i,k).lt.60000.)) then
2647                    ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2648     s              (60000.-pplay(i,k))/(60000.-30000.)
2649               else
2650                    ratqss(i,k)=ratqshaut
2651               endif
2652            enddo
2653         enddo
2654      endif
2655
2656
2657
2658
2659c  ratqs final
2660c  -----------
2661
2662      if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2663     s    .or.iflag_cldcon.ge.4) then
2664
2665! On ajoute une constante au ratqsc*2 pour tenir compte de
2666! fluctuations turbulentes de petite echelle
2667
2668         do k=1,klev
2669            do i=1,klon
2670               if ((fm_therm(i,k).gt.1.e-10)) then
2671                  ratqsc(i,k)=sqrt(ratqsc(i,k)**2+0.05**2)
2672               endif
2673            enddo
2674         enddo
2675
2676!   les ratqs sont une combinaison de ratqss et ratqsc
2677       if(prt_level.ge.9)
2678     $       write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
2679
2680         if (tau_ratqs>1.e-10) then
2681            facteur=exp(-pdtphys/tau_ratqs)
2682         else
2683            facteur=0.
2684         endif
2685         ratqs(:,:)=ratqsc(:,:)*(1.-facteur)+ratqs(:,:)*facteur
2686!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2687! FH 22/09/2009
2688! La ligne ci-dessous faisait osciller le modele et donnait une solution
2689! assymptotique bidon et dépendant fortement du pas de temps.
2690!        ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
2691!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2692         ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
2693      else
2694!   on ne prend que le ratqs stable pour fisrtilp
2695         ratqs(:,:)=ratqss(:,:)
2696      endif
2697
2698
2699c
2700c Appeler le processus de condensation a grande echelle
2701c et le processus de precipitation
2702c-------------------------------------------------------------------------
2703      CALL fisrtilp(dtime,paprs,pplay,
2704     .           t_seri, q_seri,ptconv,ratqs,
2705     .           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2706     .           rain_lsc, snow_lsc,
2707     .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
2708     .           frac_impa, frac_nucl,
2709     .           prfl, psfl, rhcl,
2710     .           zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
2711
2712      WHERE (rain_lsc < 0) rain_lsc = 0.
2713      WHERE (snow_lsc < 0) snow_lsc = 0.
2714!-----------------------------------------------------------------------------------------
2715! ajout des tendances de la diffusion turbulente
2716      CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2717!-----------------------------------------------------------------------------------------
2718      DO k = 1, klev
2719      DO i = 1, klon
2720         cldfra(i,k) = rneb(i,k)
2721         IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2722      ENDDO
2723      ENDDO
2724      IF (check) THEN
2725         za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2726         WRITE(lunout,*)"apresilp=", za
2727         zx_t = 0.0
2728         za = 0.0
2729         DO i = 1, klon
2730            za = za + airephy(i)/REAL(klon)
2731            zx_t = zx_t + (rain_lsc(i)
2732     .                  + snow_lsc(i))*airephy(i)/REAL(klon)
2733        ENDDO
2734         zx_t = zx_t/za*dtime
2735         WRITE(lunout,*)"Precip=", zx_t
2736      ENDIF
2737cIM
2738      IF (ip_ebil_phy.ge.2) THEN
2739        ztit='after fisrt'
2740        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2741     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2742     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2743        call diagphy(airephy,ztit,ip_ebil_phy
2744     e      , zero_v, zero_v, zero_v, zero_v, zero_v
2745     e      , zero_v, rain_lsc, snow_lsc, ztsol
2746     e      , d_h_vcol, d_qt, d_ec
2747     s      , fs_bound, fq_bound )
2748      END IF
2749
2750      if (mydebug) then
2751        call writefield_phy('u_seri',u_seri,llm)
2752        call writefield_phy('v_seri',v_seri,llm)
2753        call writefield_phy('t_seri',t_seri,llm)
2754        call writefield_phy('q_seri',q_seri,llm)
2755      endif
2756
2757c
2758c-------------------------------------------------------------------
2759c  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2760c-------------------------------------------------------------------
2761
2762c 1. NUAGES CONVECTIFS
2763c
2764cIM cf FH
2765c     IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2766      IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2767       snow_tiedtke=0.
2768c     print*,'avant calcul de la pseudo precip '
2769c     print*,'iflag_cldcon',iflag_cldcon
2770       if (iflag_cldcon.eq.-1) then
2771          rain_tiedtke=rain_con
2772       else
2773c       print*,'calcul de la pseudo precip '
2774          rain_tiedtke=0.
2775c         print*,'calcul de la pseudo precip 0'
2776          do k=1,klev
2777          do i=1,klon
2778             if (d_q_con(i,k).lt.0.) then
2779                rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2780     s         *(paprs(i,k)-paprs(i,k+1))/rg
2781             endif
2782          enddo
2783          enddo
2784       endif
2785c
2786c     call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2787c
2788
2789c Nuages diagnostiques pour Tiedtke
2790      CALL diagcld1(paprs,pplay,
2791cIM cf FH  .             rain_con,snow_con,ibas_con,itop_con,
2792     .             rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2793     .             diafra,dialiq)
2794      DO k = 1, klev
2795      DO i = 1, klon
2796      IF (diafra(i,k).GT.cldfra(i,k)) THEN
2797         cldliq(i,k) = dialiq(i,k)
2798         cldfra(i,k) = diafra(i,k)
2799      ENDIF
2800      ENDDO
2801      ENDDO
2802
2803      ELSE IF (iflag_cldcon.ge.3) THEN
2804c  On prend pour les nuages convectifs le max du calcul de la
2805c  convection et du calcul du pas de temps precedent diminue d'un facteur
2806c  facttemps
2807      facteur = pdtphys *facttemps
2808      do k=1,klev
2809         do i=1,klon
2810            rnebcon(i,k)=rnebcon(i,k)*facteur
2811            if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k))
2812     s      then
2813                rnebcon(i,k)=rnebcon0(i,k)
2814                clwcon(i,k)=clwcon0(i,k)
2815            endif
2816         enddo
2817      enddo
2818
2819c
2820cjq - introduce the aerosol direct and first indirect radiative forcings
2821cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2822      IF (ok_ade.OR.ok_aie) THEN
2823         IF (.NOT. aerosol_couple)
2824     &        CALL readaerosol_optic(
2825     &        debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref,
2826     &        pdtphys, pplay, paprs, t_seri, rhcl, presnivs,
2827     &        mass_solu_aero, mass_solu_aero_pi,
2828     &        tau_aero, piz_aero, cg_aero,
2829     &        tausum_aero, tau3d_aero)
2830      ELSE
2831cIM 170310 BEG
2832         tausum_aero(:,:,:) = 0.
2833cIM 170310 END
2834         tau_aero(:,:,:,:) = 0.
2835         piz_aero(:,:,:,:) = 0.
2836         cg_aero(:,:,:,:)  = 0.
2837      ENDIF
2838
2839cIM calcul nuages par le simulateur ISCCP
2840c
2841#ifdef histISCCP
2842      IF (ok_isccp) THEN
2843c
2844cIM lecture invtau, tautab des fichiers formattes
2845c
2846      IF (debut) THEN
2847c$OMP MASTER
2848c
2849      open(99,file='tautab.formatted', FORM='FORMATTED')
2850      read(99,'(f30.20)') tautab_omp
2851      close(99)
2852c
2853      open(99,file='invtau.formatted',form='FORMATTED')
2854      read(99,'(i10)') invtau_omp
2855
2856c     print*,'calcul_simulISCCP invtau_omp',invtau_omp
2857c     write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
2858
2859      close(99)
2860c$OMP END MASTER
2861c$OMP BARRIER
2862      tautab=tautab_omp
2863      invtau=invtau_omp
2864c
2865      ENDIF !debut
2866c
2867cIM appel simulateur toutes les  NINT(freq_ISCCP/dtime) heures
2868       IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2869#include "calcul_simulISCCP.h"
2870       ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2871      ENDIF !ok_isccp
2872#endif
2873
2874c   On prend la somme des fractions nuageuses et des contenus en eau
2875      cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2876      cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2877
2878      ENDIF
2879c
2880c 2. NUAGES STARTIFORMES
2881c
2882      IF (ok_stratus) THEN
2883      CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2884      DO k = 1, klev
2885      DO i = 1, klon
2886      IF (diafra(i,k).GT.cldfra(i,k)) THEN
2887         cldliq(i,k) = dialiq(i,k)
2888         cldfra(i,k) = diafra(i,k)
2889      ENDIF
2890      ENDDO
2891      ENDDO
2892      ENDIF
2893c
2894c Precipitation totale
2895c
2896      DO i = 1, klon
2897         rain_fall(i) = rain_con(i) + rain_lsc(i)
2898         snow_fall(i) = snow_con(i) + snow_lsc(i)
2899      ENDDO
2900cIM
2901      IF (ip_ebil_phy.ge.2) THEN
2902        ztit="after diagcld"
2903        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2904     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2905     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2906        call diagphy(airephy,ztit,ip_ebil_phy
2907     e      , zero_v, zero_v, zero_v, zero_v, zero_v
2908     e      , zero_v, zero_v, zero_v, ztsol
2909     e      , d_h_vcol, d_qt, d_ec
2910     s      , fs_bound, fq_bound )
2911      END IF
2912c
2913c Calculer l'humidite relative pour diagnostique
2914c
2915      DO k = 1, klev
2916      DO i = 1, klon
2917         zx_t = t_seri(i,k)
2918         IF (thermcep) THEN
2919            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2920            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2921            zx_qs  = MIN(0.5,zx_qs)
2922            zcor   = 1./(1.-retv*zx_qs)
2923            zx_qs  = zx_qs*zcor
2924         ELSE
2925           IF (zx_t.LT.t_coup) THEN
2926              zx_qs = qsats(zx_t)/pplay(i,k)
2927           ELSE
2928              zx_qs = qsatl(zx_t)/pplay(i,k)
2929           ENDIF
2930         ENDIF
2931         zx_rh(i,k) = q_seri(i,k)/zx_qs
2932         zqsat(i,k)=zx_qs
2933      ENDDO
2934      ENDDO
2935
2936cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2937c   equivalente a 2m (tpote) pour diagnostique
2938c
2939      DO i = 1, klon
2940       tpot(i)=zt2m(i)*(100000./paprs(i,1))**RKAPPA
2941       IF (thermcep) THEN
2942        IF(zt2m(i).LT.RTT) then
2943         Lheat=RLSTT
2944        ELSE
2945         Lheat=RLVTT
2946        ENDIF
2947       ELSE
2948        IF (zt2m(i).LT.RTT) THEN
2949         Lheat=RLSTT
2950        ELSE
2951         Lheat=RLVTT
2952        ENDIF
2953       ENDIF
2954       tpote(i) = tpot(i)*     
2955     . EXP((Lheat *qsat2m(i))/(RCPD*zt2m(i)))
2956      ENDDO
2957
2958      IF (config_inca /= 'none') THEN
2959#ifdef INCA
2960         CALL VTe(VTphysiq)
2961         CALL VTb(VTinca)
2962         calday = REAL(days_elapsed + 1) + jH_cur
2963
2964         call chemtime(itap+itau_phy-1, date0, dtime)
2965         IF (config_inca == 'aero') THEN
2966            CALL AEROSOL_METEO_CALC(
2967     $           calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2968     $           prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
2969         END IF
2970
2971         zxsnow_dummy(:) = 0.0
2972
2973         CALL chemhook_begin (calday,
2974     $                          days_elapsed+1,
2975     $                          jH_cur,
2976     $                          pctsrf(1,1),
2977     $                          rlat,
2978     $                          rlon,
2979     $                          airephy,
2980     $                          paprs,
2981     $                          pplay,
2982     $                          coefh,
2983     $                          pphi,
2984     $                          t_seri,
2985     $                          u,
2986     $                          v,
2987     $                          wo(:, :, 1),
2988     $                          q_seri,
2989     $                          zxtsol,
2990     $                          zxsnow_dummy,
2991     $                          solsw,
2992     $                          albsol1,
2993     $                          rain_fall,
2994     $                          snow_fall,
2995     $                          itop_con,
2996     $                          ibas_con,
2997     $                          cldfra,
2998     $                          iim,
2999     $                          jjm,
3000     $                          tr_seri,
3001     $                          ftsol,
3002     $                          paprs,
3003     $                          cdragh,
3004     $                          cdragm,
3005     $                          pctsrf,
3006     $                          pdtphys,
3007     $                            itap)
3008
3009         CALL VTe(VTinca)
3010         CALL VTb(VTphysiq)
3011#endif
3012      END IF !config_inca /= 'none'
3013c     
3014c Calculer les parametres optiques des nuages et quelques
3015c parametres pour diagnostiques:
3016c
3017
3018      IF (aerosol_couple) THEN
3019         mass_solu_aero(:,:)    = ccm(:,:,1)
3020         mass_solu_aero_pi(:,:) = ccm(:,:,2)
3021      END IF
3022
3023      if (ok_newmicro) then
3024      CALL newmicro (paprs, pplay,ok_newmicro,
3025     .            t_seri, cldliq, cldfra, cldtau, cldemi,
3026     .            cldh, cldl, cldm, cldt, cldq,
3027     .            flwp, fiwp, flwc, fiwc,
3028     e            ok_aie,
3029     e            mass_solu_aero, mass_solu_aero_pi,
3030     e            bl95_b0, bl95_b1,
3031     s            cldtaupi, re, fl, ref_liq, ref_ice)
3032      else
3033      CALL nuage (paprs, pplay,
3034     .            t_seri, cldliq, cldfra, cldtau, cldemi,
3035     .            cldh, cldl, cldm, cldt, cldq,
3036     e            ok_aie,
3037     e            mass_solu_aero, mass_solu_aero_pi,
3038     e            bl95_b0, bl95_b1,
3039     s            cldtaupi, re, fl)
3040     
3041      endif
3042c
3043c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3044c
3045      IF (MOD(itaprad,radpas).EQ.0) THEN
3046
3047      DO i = 1, klon
3048         albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
3049     .             + falb1(i,is_lic) * pctsrf(i,is_lic)
3050     .             + falb1(i,is_ter) * pctsrf(i,is_ter)
3051     .             + falb1(i,is_sic) * pctsrf(i,is_sic)
3052         albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
3053     .               + falb2(i,is_lic) * pctsrf(i,is_lic)
3054     .               + falb2(i,is_ter) * pctsrf(i,is_ter)
3055     .               + falb2(i,is_sic) * pctsrf(i,is_sic)
3056      ENDDO
3057
3058      if (mydebug) then
3059        call writefield_phy('u_seri',u_seri,llm)
3060        call writefield_phy('v_seri',v_seri,llm)
3061        call writefield_phy('t_seri',t_seri,llm)
3062       call writefield_phy('q_seri',q_seri,llm)
3063      endif
3064     
3065      IF (aerosol_couple) THEN
3066#ifdef INCA
3067         CALL radlwsw_inca
3068     e        (kdlon,kflev,dist, rmu0, fract, solaire,
3069     e        paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
3070     e        wo(:, :, 1),
3071     e        cldfra, cldemi, cldtau,
3072     s        heat,heat0,cool,cool0,radsol,albpla,
3073     s        topsw,toplw,solsw,sollw,
3074     s        sollwdown,
3075     s        topsw0,toplw0,solsw0,sollw0,
3076     s        lwdn0, lwdn, lwup0, lwup,
3077     s        swdn0, swdn, swup0, swup,
3078     e        ok_ade, ok_aie,
3079     e        tau_aero, piz_aero, cg_aero,
3080     s        topswad_aero, solswad_aero,
3081     s        topswad0_aero, solswad0_aero,
3082     s        topsw_aero, topsw0_aero,
3083     s        solsw_aero, solsw0_aero,
3084     e        cldtaupi,
3085     s        topswai_aero, solswai_aero)
3086           
3087#endif
3088      ELSE
3089
3090         CALL radlwsw
3091     e        (dist, rmu0, fract,
3092     e        paprs, pplay,zxtsol,albsol1, albsol2,
3093     e        t_seri,q_seri,wo,
3094     e        cldfra, cldemi, cldtau,
3095     e        ok_ade, ok_aie,
3096     e        tau_aero, piz_aero, cg_aero,
3097     e        cldtaupi,new_aod,
3098     e        zqsat, flwc, fiwc,
3099     s        heat,heat0,cool,cool0,radsol,albpla,
3100     s        topsw,toplw,solsw,sollw,
3101     s        sollwdown,
3102     s        topsw0,toplw0,solsw0,sollw0,
3103     s        lwdn0, lwdn, lwup0, lwup,
3104     s        swdn0, swdn, swup0, swup,
3105     s        topswad_aero, solswad_aero,
3106     s        topswai_aero, solswai_aero,
3107     o        topswad0_aero, solswad0_aero,
3108     o        topsw_aero, topsw0_aero,
3109     o        solsw_aero, solsw0_aero,
3110     o        topswcf_aero, solswcf_aero)
3111         
3112
3113      ENDIF ! aerosol_couple
3114      itaprad = 0
3115      ENDIF ! MOD(itaprad,radpas)
3116      itaprad = itaprad + 1
3117
3118      if (iflag_radia.eq.0 .and. prt_level.ge.9) then
3119      print *,'--------------------------------------------------'
3120      print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3121      print *,'>>>>           heat et cool mis a zero '
3122      print *,'--------------------------------------------------'
3123      heat=0.
3124      cool=0.
3125      endif
3126
3127c
3128c Ajouter la tendance des rayonnements (tous les pas)
3129c
3130      DO k = 1, klev
3131      DO i = 1, klon
3132         t_seri(i,k) = t_seri(i,k)
3133     .               + (heat(i,k)-cool(i,k)) * dtime/RDAY
3134      ENDDO
3135      ENDDO
3136c
3137      if (mydebug) then
3138        call writefield_phy('u_seri',u_seri,llm)
3139        call writefield_phy('v_seri',v_seri,llm)
3140        call writefield_phy('t_seri',t_seri,llm)
3141        call writefield_phy('q_seri',q_seri,llm)
3142      endif
3143 
3144cIM
3145      IF (ip_ebil_phy.ge.2) THEN
3146        ztit='after rad'
3147        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3148     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3149     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3150        call diagphy(airephy,ztit,ip_ebil_phy
3151     e      , topsw, toplw, solsw, sollw, zero_v
3152     e      , zero_v, zero_v, zero_v, ztsol
3153     e      , d_h_vcol, d_qt, d_ec
3154     s      , fs_bound, fq_bound )
3155      END IF
3156c
3157c
3158c Calculer l'hydrologie de la surface
3159c
3160c      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3161c     .            agesno, ftsol,fqsurf,fsnow, ruis)
3162c
3163
3164c
3165c Calculer le bilan du sol et la derive de temperature (couplage)
3166c
3167      DO i = 1, klon
3168c         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3169c a la demande de JLD
3170         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3171      ENDDO
3172c
3173cmoddeblott(jan95)
3174c Appeler le programme de parametrisation de l'orographie
3175c a l'echelle sous-maille:
3176c
3177      IF (ok_orodr) THEN
3178c
3179c  selection des points pour lesquels le shema est actif:
3180        igwd=0
3181        DO i=1,klon
3182        itest(i)=0
3183c        IF ((zstd(i).gt.10.0)) THEN
3184        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3185          itest(i)=1
3186          igwd=igwd+1
3187          idx(igwd)=i
3188        ENDIF
3189        ENDDO
3190c        igwdim=MAX(1,igwd)
3191c
3192        IF (ok_strato) THEN
3193       
3194          CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
3195     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
3196     e                   igwd,idx,itest,
3197     e                   t_seri, u_seri, v_seri,
3198     s                   zulow, zvlow, zustrdr, zvstrdr,
3199     s                   d_t_oro, d_u_oro, d_v_oro)
3200
3201       ELSE
3202        CALL drag_noro(klon,klev,dtime,paprs,pplay,
3203     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
3204     e                   igwd,idx,itest,
3205     e                   t_seri, u_seri, v_seri,
3206     s                   zulow, zvlow, zustrdr, zvstrdr,
3207     s                   d_t_oro, d_u_oro, d_v_oro)
3208       ENDIF
3209c
3210c  ajout des tendances
3211!-----------------------------------------------------------------------------------------
3212! ajout des tendances de la trainee de l'orographie
3213      CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3214!-----------------------------------------------------------------------------------------
3215c
3216      ENDIF ! fin de test sur ok_orodr
3217c
3218      if (mydebug) then
3219        call writefield_phy('u_seri',u_seri,llm)
3220        call writefield_phy('v_seri',v_seri,llm)
3221        call writefield_phy('t_seri',t_seri,llm)
3222        call writefield_phy('q_seri',q_seri,llm)
3223      endif
3224     
3225      IF (ok_orolf) THEN
3226c
3227c  selection des points pour lesquels le shema est actif:
3228        igwd=0
3229        DO i=1,klon
3230        itest(i)=0
3231        IF ((zpic(i)-zmea(i)).GT.100.) THEN
3232          itest(i)=1
3233          igwd=igwd+1
3234          idx(igwd)=i
3235        ENDIF
3236        ENDDO
3237c        igwdim=MAX(1,igwd)
3238c
3239        IF (ok_strato) THEN
3240
3241          CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
3242     e                   rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
3243     e                   igwd,idx,itest,
3244     e                   t_seri, u_seri, v_seri,
3245     s                   zulow, zvlow, zustrli, zvstrli,
3246     s                   d_t_lif, d_u_lif, d_v_lif               )
3247       
3248        ELSE
3249          CALL lift_noro(klon,klev,dtime,paprs,pplay,
3250     e                   rlat,zmea,zstd,zpic,
3251     e                   itest,
3252     e                   t_seri, u_seri, v_seri,
3253     s                   zulow, zvlow, zustrli, zvstrli,
3254     s                   d_t_lif, d_u_lif, d_v_lif)
3255       ENDIF
3256c   
3257!-----------------------------------------------------------------------------------------
3258! ajout des tendances de la portance de l'orographie
3259      CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3260!-----------------------------------------------------------------------------------------
3261c
3262      ENDIF ! fin de test sur ok_orolf
3263C  HINES GWD PARAMETRIZATION
3264
3265       IF (ok_hines) then
3266
3267         CALL hines_gwd(klon,klev,dtime,paprs,pplay,
3268     i                  rlat,t_seri,u_seri,v_seri,
3269     o                  zustrhi,zvstrhi,
3270     o                  d_t_hin, d_u_hin, d_v_hin)
3271c
3272c  ajout des tendances
3273        CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'lif')
3274
3275      ENDIF
3276c
3277
3278c
3279cIM cf. FLott BEG
3280C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3281
3282      if (mydebug) then
3283        call writefield_phy('u_seri',u_seri,llm)
3284        call writefield_phy('v_seri',v_seri,llm)
3285        call writefield_phy('t_seri',t_seri,llm)
3286        call writefield_phy('q_seri',q_seri,llm)
3287      endif
3288
3289      DO i = 1, klon
3290        zustrph(i)=0.
3291        zvstrph(i)=0.
3292      ENDDO
3293      DO k = 1, klev
3294      DO i = 1, klon
3295       zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3296     c            (paprs(i,k)-paprs(i,k+1))/rg
3297       zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3298     c            (paprs(i,k)-paprs(i,k+1))/rg
3299      ENDDO
3300      ENDDO
3301c
3302cIM calcul composantes axiales du moment angulaire et couple des montagnes
3303c
3304      IF (is_sequential) THEN
3305     
3306        CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
3307     C                 ra,rg,romega,
3308     C                 rlat,rlon,pphis,
3309     C                 zustrdr,zustrli,zustrph,
3310     C                 zvstrdr,zvstrli,zvstrph,
3311     C                 paprs,u,v,
3312     C                 aam, torsfc)
3313       ENDIF
3314cIM cf. FLott END
3315cIM
3316      IF (ip_ebil_phy.ge.2) THEN
3317        ztit='after orography'
3318        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3319     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3320     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3321         call diagphy(airephy,ztit,ip_ebil_phy
3322     e      , zero_v, zero_v, zero_v, zero_v, zero_v
3323     e      , zero_v, zero_v, zero_v, ztsol
3324     e      , d_h_vcol, d_qt, d_ec
3325     s      , fs_bound, fq_bound )
3326      END IF
3327c
3328c
3329!====================================================================
3330! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3331!====================================================================
3332! Abderrahmane 24.08.09
3333
3334      IF (ok_cosp) THEN
3335! adeclarer
3336#ifdef CPP_COSP
3337       IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3338
3339       print*,'freq_cosp',freq_cosp
3340          mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3341!       print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3342!     s        ref_liq,ref_ice
3343          call phys_cosp(itap,dtime,freq_cosp,
3344     $                 ecrit_mth,ecrit_day,ecrit_hf,overlap,
3345     $                   klon,klev,rlon,rlat,presnivs,
3346     $                   ref_liq,ref_ice,
3347     $                   pctsrf(:,is_ter)+pctsrf(:,is_lic),
3348     $                   zu10m,zv10m,
3349     $                   zphi,paprs(:,1:klev),pplay,zxtsol,t_seri,
3350     $                   qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc,
3351     $                   prfl(:,1:klev),psfl(:,1:klev),
3352     $                   pmflxr(:,1:klev),pmflxs(:,1:klev),
3353     $                   mr_ozone,cldtau, cldemi)
3354!     L          calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3355!     L          cfaddbze,clcalipso2,dbze,cltlidarradar,
3356!     M          clMISR,
3357!     R          clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3358!     I          tauisccp,albisccp,meantbisccp,meantbclrisccp)
3359
3360         ENDIF
3361
3362#endif
3363       ENDIF  !ok_cosp
3364!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3365cAA
3366cAA Installation de l'interface online-offline pour traceurs
3367cAA
3368c====================================================================
3369c   Calcul  des tendances traceurs
3370c====================================================================
3371C
3372
3373      call phytrac (
3374     I     itap,     days_elapsed+1,    jH_cur,   debut,
3375     I     lafin,    dtime,     u, v,     t,
3376     I     paprs,    pplay,     pmfu,     pmfd,
3377     I     pen_u,    pde_u,     pen_d,    pde_d,
3378     I     cdragh,   coefh,     fm_therm, entr_therm,
3379     I     u1,       v1,        ftsol,    pctsrf,
3380     I     rlat,     frac_impa, frac_nucl,rlon,
3381     I     presnivs, pphis,     pphi,     albsol1,
3382     I     qx(:,:,ivap),rhcl,   cldfra,   rneb,
3383     I     diafra,   cldliq,    itop_con, ibas_con,
3384     I     pmflxr,   pmflxs,    prfl,     psfl,
3385     I     da,       phi,       mp,       upwd,     
3386     I     dnwd,     aerosol_couple,      flxmass_w,
3387     I     tau_aero, piz_aero,  cg_aero,  ccm,
3388     I     rfname,
3389     O     tr_seri)
3390
3391      IF (offline) THEN
3392
3393       IF (prt_level.ge.9)
3394     $    print*,'Attention on met a 0 les thermiques pour phystoke'
3395         call phystokenc (
3396     I                   nlon,klev,pdtphys,rlon,rlat,
3397     I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3398     I                   fm_therm,entr_therm,
3399     I                   cdragh,coefh,u1,v1,ftsol,pctsrf,
3400     I                   frac_impa, frac_nucl,
3401     I                   pphis,airephy,dtime,itap)
3402
3403
3404      ENDIF
3405
3406c
3407c Calculer le transport de l'eau et de l'energie (diagnostique)
3408c
3409      CALL transp (paprs,zxtsol,
3410     e                   t_seri, q_seri, u_seri, v_seri, zphi,
3411     s                   ve, vq, ue, uq)
3412c
3413cIM global posePB BEG
3414      IF(1.EQ.0) THEN
3415c
3416      CALL transp_lay (paprs,zxtsol,
3417     e                   t_seri, q_seri, u_seri, v_seri, zphi,
3418     s                   ve_lay, vq_lay, ue_lay, uq_lay)
3419c
3420      ENDIF !(1.EQ.0) THEN
3421cIM global posePB END
3422c Accumuler les variables a stocker dans les fichiers histoire:
3423c
3424c+jld ec_conser
3425      DO k = 1, klev
3426      DO i = 1, klon
3427        ZRCPD = RCPD*(1.0+RVTMP2*q_seri(i,k))
3428        d_t_ec(i,k)=0.5/ZRCPD
3429     $      *(u(i,k)**2+v(i,k)**2-u_seri(i,k)**2-v_seri(i,k)**2)
3430      ENDDO
3431      ENDDO
3432
3433      DO k = 1, klev
3434      DO i = 1, klon
3435        t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3436        d_t_ec(i,k) = d_t_ec(i,k)/dtime
3437       END DO
3438      END DO
3439c-jld ec_conser
3440cIM
3441      IF (ip_ebil_phy.ge.1) THEN
3442        ztit='after physic'
3443        CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3444     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3445     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3446C     Comme les tendances de la physique sont ajoute dans la dynamique,
3447C     on devrait avoir que la variation d'entalpie par la dynamique
3448C     est egale a la variation de la physique au pas de temps precedent.
3449C     Donc la somme de ces 2 variations devrait etre nulle.
3450
3451        call diagphy(airephy,ztit,ip_ebil_phy
3452     e      , topsw, toplw, solsw, sollw, sens
3453     e      , evap, rain_fall, snow_fall, ztsol
3454     e      , d_h_vcol, d_qt, d_ec
3455     s      , fs_bound, fq_bound )
3456C
3457      d_h_vcol_phy=d_h_vcol
3458C
3459      END IF
3460C
3461c=======================================================================
3462c   SORTIES
3463c=======================================================================
3464
3465cIM Interpolation sur les niveaux de pression du NMC
3466c   -------------------------------------------------
3467c
3468#include "calcul_STDlev.h"
3469      twriteSTD(:,:,1)=tsumSTD(:,:,1)
3470      qwriteSTD(:,:,1)=qsumSTD(:,:,1)
3471      rhwriteSTD(:,:,1)=rhsumSTD(:,:,1)
3472      phiwriteSTD(:,:,1)=phisumSTD(:,:,1)
3473      uwriteSTD(:,:,1)=usumSTD(:,:,1)
3474      vwriteSTD(:,:,1)=vsumSTD(:,:,1)
3475      wwriteSTD(:,:,1)=wsumSTD(:,:,1)
3476
3477      twriteSTD(:,:,2)=tsumSTD(:,:,2)
3478      qwriteSTD(:,:,2)=qsumSTD(:,:,2)
3479      rhwriteSTD(:,:,2)=rhsumSTD(:,:,2)
3480      phiwriteSTD(:,:,2)=phisumSTD(:,:,2)
3481      uwriteSTD(:,:,2)=usumSTD(:,:,2)
3482      vwriteSTD(:,:,2)=vsumSTD(:,:,2)
3483      wwriteSTD(:,:,2)=wsumSTD(:,:,2)
3484
3485      twriteSTD(:,:,3)=tlevSTD(:,:)
3486      qwriteSTD(:,:,3)=qlevSTD(:,:)
3487      rhwriteSTD(:,:,3)=rhlevSTD(:,:)
3488      phiwriteSTD(:,:,3)=philevSTD(:,:)
3489      uwriteSTD(:,:,3)=ulevSTD(:,:)
3490      vwriteSTD(:,:,3)=vlevSTD(:,:)
3491      wwriteSTD(:,:,3)=wlevSTD(:,:)
3492
3493      twriteSTD(:,:,4)=tlevSTD(:,:)
3494      qwriteSTD(:,:,4)=qlevSTD(:,:)
3495      rhwriteSTD(:,:,4)=rhlevSTD(:,:)
3496      phiwriteSTD(:,:,4)=philevSTD(:,:)
3497      uwriteSTD(:,:,4)=ulevSTD(:,:)
3498      vwriteSTD(:,:,4)=vlevSTD(:,:)
3499      wwriteSTD(:,:,4)=wlevSTD(:,:)
3500c
3501cIM ajoute 5eme niveau 170310 BEG
3502      twriteSTD(:,:,5)=tlevSTD(:,:)
3503      qwriteSTD(:,:,5)=qlevSTD(:,:)
3504      rhwriteSTD(:,:,5)=rhlevSTD(:,:)
3505      phiwriteSTD(:,:,5)=philevSTD(:,:)
3506      uwriteSTD(:,:,5)=ulevSTD(:,:)
3507      vwriteSTD(:,:,5)=vlevSTD(:,:)
3508      wwriteSTD(:,:,5)=wlevSTD(:,:)
3509cIM for NMC files
3510      DO n=1, nlevSTD3
3511       DO k=1, nlevSTD
3512        if(rlevSTD3(n).EQ.rlevSTD(k)) THEN
3513         twriteSTD3(:,n)=tlevSTD(:,k)
3514         qwriteSTD3(:,n)=qlevSTD(:,k)
3515         rhwriteSTD3(:,n)=rhlevSTD(:,k)
3516         phiwriteSTD3(:,n)=philevSTD(:,k)
3517         uwriteSTD3(:,n)=ulevSTD(:,k)
3518         vwriteSTD3(:,n)=vlevSTD(:,k)
3519         wwriteSTD3(:,n)=wlevSTD(:,k)
3520        endif !rlevSTD3(n).EQ.rlevSTD(k)
3521       ENDDO
3522      ENDDO
3523c
3524      DO n=1, nlevSTD8
3525       DO k=1, nlevSTD
3526        if(rlevSTD8(n).EQ.rlevSTD(k)) THEN
3527         tnondefSTD8(:,n)=tnondef(:,k,2)
3528         twriteSTD8(:,n)=tsumSTD(:,k,2)
3529         qwriteSTD8(:,n)=qsumSTD(:,k,2)
3530         rhwriteSTD8(:,n)=rhsumSTD(:,k,2)
3531         phiwriteSTD8(:,n)=phisumSTD(:,k,2)
3532         uwriteSTD8(:,n)=usumSTD(:,k,2)
3533         vwriteSTD8(:,n)=vsumSTD(:,k,2)
3534         wwriteSTD8(:,n)=wsumSTD(:,k,2)
3535        endif !rlevSTD8(n).EQ.rlevSTD(k)
3536       ENDDO
3537      ENDDO
3538c
3539c slp sea level pressure
3540      slp(:) = paprs(:,1)*exp(pphis(:)/(RD*t_seri(:,1)))
3541c
3542ccc prw = eau precipitable
3543      DO i = 1, klon
3544       prw(i) = 0.
3545       DO k = 1, klev
3546        prw(i) = prw(i) +
3547     .           q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3548       ENDDO
3549      ENDDO
3550c
3551cIM initialisation + calculs divers diag AMIP2
3552c
3553#include "calcul_divers.h"
3554c
3555      IF (config_inca /= 'none') THEN
3556#ifdef INCA
3557         CALL VTe(VTphysiq)
3558         CALL VTb(VTinca)
3559
3560         CALL chemhook_end (
3561     $                        dtime,
3562     $                        pplay,
3563     $                        t_seri,
3564     $                        tr_seri,
3565     $                        nbtr,
3566     $                        paprs,
3567     $                        q_seri,
3568     $                        airephy,
3569     $                        pphi,
3570     $                        pphis,
3571     $                        zx_rh)
3572
3573         CALL VTe(VTinca)
3574         CALL VTb(VTphysiq)
3575#endif
3576      END IF
3577
3578c=============================================================
3579c
3580c Convertir les incrementations en tendances
3581c
3582      if (mydebug) then
3583        call writefield_phy('u_seri',u_seri,llm)
3584        call writefield_phy('v_seri',v_seri,llm)
3585        call writefield_phy('t_seri',t_seri,llm)
3586        call writefield_phy('q_seri',q_seri,llm)
3587      endif
3588
3589      DO k = 1, klev
3590      DO i = 1, klon
3591         d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3592         d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3593         d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3594         d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3595         d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3596      ENDDO
3597      ENDDO
3598c
3599      IF (nqtot.GE.3) THEN
3600      DO iq = 3, nqtot
3601      DO  k = 1, klev
3602      DO  i = 1, klon
3603         d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3604      ENDDO
3605      ENDDO
3606      ENDDO
3607      ENDIF
3608c
3609cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3610cIM global posePB#include "write_bilKP_ins.h"
3611cIM global posePB#include "write_bilKP_ave.h"
3612c
3613c Sauvegarder les valeurs de t et q a la fin de la physique:
3614c
3615      DO k = 1, klev
3616      DO i = 1, klon
3617         u_ancien(i,k) = u_seri(i,k)
3618         v_ancien(i,k) = v_seri(i,k)
3619         t_ancien(i,k) = t_seri(i,k)
3620         q_ancien(i,k) = q_seri(i,k)
3621      ENDDO
3622      ENDDO
3623c
3624!==========================================================================
3625! Sorties des tendances pour un point particulier
3626! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3627! pour le debug
3628! La valeur de igout est attribuee plus haut dans le programme
3629!==========================================================================
3630
3631      if (prt_level.ge.1) then
3632      write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3633      write(lunout,*)
3634     s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
3635      write(lunout,*)
3636     s  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
3637     s  pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
3638     s  pctsrf(igout,is_sic)
3639      write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3640      do k=1,klev
3641         write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
3642     s   d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
3643     s   d_t_eva(igout,k)
3644      enddo
3645      write(lunout,*) 'cool,heat'
3646      do k=1,klev
3647         write(lunout,*) cool(igout,k),heat(igout,k)
3648      enddo
3649
3650      write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3651      do k=1,klev
3652         write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
3653     s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3654      enddo
3655
3656      write(lunout,*) 'd_ps ',d_ps(igout)
3657      write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3658      do k=1,klev
3659         write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
3660     s  d_qx(igout,k,1),d_qx(igout,k,2)
3661      enddo
3662      endif
3663
3664!==========================================================================
3665
3666c============================================================
3667c   Calcul de la temperature potentielle
3668c============================================================
3669      DO k = 1, klev
3670      DO i = 1, klon
3671        theta(i,k)=t(i,k)*(100000./pplay(i,k))**(RD/RCPD)
3672      ENDDO
3673      ENDDO
3674c
3675
3676c 22.03.04 BEG
3677c=============================================================
3678c   Ecriture des sorties
3679c=============================================================
3680#ifdef CPP_IOIPSL
3681 
3682c Recupere des varibles calcule dans differents modules
3683c pour ecriture dans histxxx.nc
3684
3685      ! Get some variables from module fonte_neige_mod
3686      CALL fonte_neige_get_vars(pctsrf,
3687     .     zxfqcalving, zxfqfonte, zxffonte)
3688
3689
3690#include "phys_output_write.h"
3691
3692#ifdef histISCCP
3693#include "write_histISCCP.h"
3694#endif
3695
3696#ifdef histNMC
3697#include "write_histhfNMC.h"
3698#include "write_histdayNMC.h"
3699#include "write_histmthNMC.h"
3700#endif
3701
3702#include "write_histday_seri.h"
3703
3704#include "write_paramLMDZ_phy.h"
3705
3706#endif
3707
3708c 22.03.04 END
3709c
3710c====================================================================
3711c Si c'est la fin, il faut conserver l'etat de redemarrage
3712c====================================================================
3713c
3714     
3715
3716      IF (lafin) THEN
3717         itau_phy = itau_phy + itap
3718         CALL phyredem ("restartphy.nc")
3719!         open(97,form="unformatted",file="finbin")
3720!         write(97) u_seri,v_seri,t_seri,q_seri
3721!         close(97)
3722C$OMP MASTER
3723         if (read_climoz >= 1) then
3724            if (is_mpi_root) then
3725               call nf95_close(ncid_climoz)
3726            end if
3727            deallocate(press_climoz) ! pointer
3728         end if
3729C$OMP END MASTER
3730      ENDIF
3731     
3732!      first=.false.
3733
3734      RETURN
3735      END
3736      FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3737      IMPLICIT none
3738c
3739c Calculer et imprimer l'eau totale. A utiliser pour verifier
3740c la conservation de l'eau
3741c
3742#include "YOMCST.h"
3743      INTEGER klon,klev
3744      REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3745      REAL aire(klon)
3746      REAL qtotal, zx, qcheck
3747      INTEGER i, k
3748c
3749      zx = 0.0
3750      DO i = 1, klon
3751         zx = zx + aire(i)
3752      ENDDO
3753      qtotal = 0.0
3754      DO k = 1, klev
3755      DO i = 1, klon
3756         qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3757     .                     *(paprs(i,k)-paprs(i,k+1))/RG
3758      ENDDO
3759      ENDDO
3760c
3761      qcheck = qtotal/zx
3762c
3763      RETURN
3764      END
3765      SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3766      IMPLICIT none
3767c
3768c Tranformer une variable de la grille physique a
3769c la grille d'ecriture
3770c
3771      INTEGER nfield,nlon,iim,jjmp1, jjm
3772      REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3773c
3774      INTEGER i, n, ig
3775c
3776      jjm = jjmp1 - 1
3777      DO n = 1, nfield
3778         DO i=1,iim
3779            ecrit(i,n) = fi(1,n)
3780            ecrit(i+jjm*iim,n) = fi(nlon,n)
3781         ENDDO
3782         DO ig = 1, nlon - 2
3783           ecrit(iim+ig,n) = fi(1+ig,n)
3784         ENDDO
3785      ENDDO
3786      RETURN
3787      END
Note: See TracBrowser for help on using the repository browser.