Context Navigation

physiq.F @ 2996

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

Necessary modifications to use the model in an aquaplanet or terraplanet configuration

Modifications nécessaires pour l'utilisation du moèle en configuration aquaplanète
ou terraplanète

Dans la dynamique:
Changement pour le paramètre iflag_phys:

auparavant : 0 pas de physique, 1 phyique, 2 rappel
ici on ajoute iflag_phys>=100 pour les aqua et terra

Du coup, dans leapfrog.F et gcm.F on appelle la physique si
iflag_phys=1.or.iflag_phys>=100
Dans iniacademic, on initialise si iflag_phys>=2 au lieu de =2
Dans gcm.F, on appelle en plus iniaqua (sous une clef CPP_EARTH)
Dans iniacademic, on met ps=108080 pour les aqua et terra pour répondre
à une specification CFMIP.

Dans la physique:
On ajoute phyaqua.F qui contient :
iniaqua : initialise les startphy.nc et limit.nc pour la phyique
zenang_an : calcule un ensolleillement moyen sur l'année en fonction de

la latiude (A. Campoy).

profil_sst : qui calcule différents profils latitudinaux de SSTs.
writelim : écrit le fichier limit.nc

Dans physiq.F
On ajoute la possibilité d'appeller un calcul d'ensoleillement moyen
sur l'année quand solarlong0=1000.

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

File size: 129.7 KB

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

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

Download in other formats:

Original Format