Context Navigation

physiq.F @ 1486

Last change on this file since 1486 was 1486, checked in by Laurent Fairhead, 13 years ago

Inclusion de la routine "stats" du modele martien qui permet de sortir
le cycle diurne moyen de différentes variables et l'écart-type.
Pour activer cette routine mettre
callstats=y
dans config.def. Les résultats seront placés dans le fichier stats.nc
Pour rajouter des variables: allez à la fin de physiq.F et prendre modèle sur
les lignes call wstats(...) et ne pas oublier de modifier le nombre de variables
sorties dans statto.h (n2dvar et n2dvar).
Attention: la routine n'est pas optimale (elle stocke ses résultats intermédiaires
dans le fichier stats.nc et les lit et écrit à chaque appel de la physique)
et n'a été testée dans LMDZ pour l'instant que sur PC/gfortran
et sur SX8/monoprocesseur. Encore du travail donc. Mais la routine est utilisée
régulièrement avec le modèle martien :-)

Inclusion of the martian model "stats" routine that outputs the mean diurnal
cycle of variables as well as their standard deviation.
To activate the routine, one needs to put
callstats=y
in the config.def file. Results will be output in a stats.nc file
To add more output variables: go to the end of physiq.F and use the call wstats(...)
lines as a template and do not forget to modify the number of output variables in
statto.h (n2dvar and n3dvar)
Caution: the routine is not optimal (as it stocks preliminary results in the stats.nc
file reading and writing at each physics step) and has only been tested in LMDZ
with a PC/gfortran setup and on a NEC-SX8/monoprocessor setup. Some work then still to
be done. But it is used regularly with the martial model :-)

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 125.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format