Context Navigation

physiq.F @ 1221

Last change on this file since 1221 was 1221, checked in by jghattas, 15 years ago

Modification pour les aerosols sels marin.

Nicolas Yan, Yves Balkanski LSCE

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

File size: 114.1 KB

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

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

Download in other formats:

Original Format