Context Navigation

physiq.F @ 1669

Last change on this file since 1669 was 1669, checked in by Laurent Fairhead, 12 years ago

Version testing basée sur la r1668

http://lmdz.lmd.jussieu.fr/utilisateurs/distribution-du-modele/versions-intermediaires

Testing release based on r1668

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

File size: 141.9 KB

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

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

Download in other formats:

Original Format