source: LMDZ4/trunk/libf/phylmd/physiq.F @ 1340

Last change on this file since 1340 was 1340, checked in by Laurent Fairhead, 14 years ago

Use generic rather than specific names for intrinsic procedures.
Modifications in r1334 (on Vprecip) were not propagated to physiq.F leading to a non
running program


Utilisations des noms génériques pour les procédures intrinsèques
Les modifications sur la déclaration de Vprecip dans la révision r1334 n'ont pas été
"remontées" à physiq.F d'où problème à l'exécution (en particulier sous g95)

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