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

Last change on this file since 1422 was 1422, checked in by jghattas, 14 years ago

Corrected error in use of flag prt_level.
MPL

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