Context Navigation

physiq.F @ 969

Last change on this file since 969 was 968, checked in by Laurent Fairhead, 17 years ago
Menage un peu trop rapide sur les parametres ocean et ok_veget LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 109.3 KB

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

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

Download in other formats:

Original Format