source: LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F @ 1123

Last change on this file since 1123 was 1114, checked in by jghattas, 16 years ago

Creation du module infotrac:

  • contient les variables de advtrac.h
  • contient la subroutine iniadvtrac renommer en infotrac_init
  • le nombre des traceurs est lu dans tracer.def en dynamique (ou par default ou recu par INCA)
  • ce module est utilise dans la dynamique et la physique
  • contient aussi la variable nbtr qui avant etait stockee dans dimphy

Le fichier advtrac.h n'existe plus.
La compilation ne prend plus en compte le nombre de traceur.

/JG

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