Context Navigation

physiq.F @ 436

Last change on this file since 436 was 436, checked in by lmdzadmin, 21 years ago
Meme erreur de synchro dans l'appel a phyredem LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 63.8 KB

Line
1	c
2	c $Header$
3	c
4	SUBROUTINE physiq (nlon,nlev,nqmax,
5	. debut,lafin,rjourvrai,rjour_ecri,gmtime,pdtphys,
6	. paprs,pplay,pphi,pphis,paire,presnivs,clesphy0,
7	. u,v,t,qx,
8	. omega, cufi, cvfi,
9	. d_u, d_v, d_t, d_qx, d_ps)
10	USE ioipsl
11	USE histcom
12	USE writephys
13
14	IMPLICIT none
15	c======================================================================
16	c
17	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
18	c
19	c Objet: Moniteur general de la physique du modele
20	cAA Modifications quant aux traceurs :
21	cAA - uniformisation des parametrisations ds phytrac
22	cAA - stockage des moyennes des champs necessaires
23	cAA en mode traceur off-line
24	c======================================================================
25	c modif ( P. Le Van , 12/10/98 )
26	c
27	c Arguments:
28	c
29	c nlon----input-I-nombre de points horizontaux
30	c nlev----input-I-nombre de couches verticales
31	c nqmax---input-I-nombre de traceurs (y compris vapeur d'eau) = 1
32	c debut---input-L-variable logique indiquant le premier passage
33	c lafin---input-L-variable logique indiquant le dernier passage
34	c rjour---input-R-numero du jour de l'experience
35	c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
36	c pdtphys-input-R-pas d'integration pour la physique (seconde)
37	c paprs---input-R-pression pour chaque inter-couche (en Pa)
38	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
39	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
40	c pphis---input-R-geopotentiel du sol
41	c paire---input-R-aire de chaque maille
42	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
43	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
44	c v-------input-R-vitesse Y (de S a N) en m/s
45	c t-------input-R-temperature (K)
46	c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
47	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
48	c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
49	c omega---input-R-vitesse verticale en Pa/s
50	c cufi----input-R-resolution des mailles en x (m)
51	c cvfi----input-R-resolution des mailles en y (m)
52	c
53	c d_u-----output-R-tendance physique de "u" (m/s/s)
54	c d_v-----output-R-tendance physique de "v" (m/s/s)
55	c d_t-----output-R-tendance physique de "t" (K/s)
56	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
57	c d_ps----output-R-tendance physique de la pression au sol
58	c======================================================================
59	#include "dimensions.h"
60	integer jjmp1
61	parameter (jjmp1=jjm+1-1/jjm)
62	#include "dimphy.h"
63	#include "regdim.h"
64	#include "indicesol.h"
65	#include "dimsoil.h"
66	#include "clesphys.h"
67	#include "control.h"
68	#include "temps.h"
69	c======================================================================
70	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
71	PARAMETER (ok_cvl=.TRUE.)
72	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
73	PARAMETER (ok_gust=.FALSE.)
74	c======================================================================
75	LOGICAL check ! Verifier la conservation du modele en eau
76	PARAMETER (check=.FALSE.)
77	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
78	PARAMETER (ok_stratus=.FALSE.)
79	c======================================================================
80	c Parametres lies au coupleur OASIS:
81	#include "oasis.h"
82	INTEGER,SAVE :: npas, nexca
83	logical rnpb
84	parameter(rnpb=.true.)
85	c PARAMETER (npas=1440)
86	c PARAMETER (nexca=48)
87	EXTERNAL fromcpl, intocpl, inicma
88	c ocean = type de modele ocean a utiliser: force, slab, couple
89	character*6 ocean
90	SAVE ocean
91
92	c parameter (ocean = 'force ')
93	c parameter (ocean = 'couple')
94	logical ok_ocean
95	c======================================================================
96	c Clef controlant l'activation du cycle diurne:
97	ccc LOGICAL cycle_diurne
98	ccc PARAMETER (cycle_diurne=.FALSE.)
99	c======================================================================
100	c Modele thermique du sol, a activer pour le cycle diurne:
101	ccc LOGICAL soil_model
102	ccc PARAMETER (soil_model=.FALSE.)
103	logical ok_veget
104	save ok_veget
105	c parameter (ok_veget = .true.)
106	c parameter (ok_veget = .false.)
107	c======================================================================
108	c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
109	c le calcul du rayonnement est celle apres la precipitation des nuages.
110	c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
111	c la condensation et la precipitation. Cette cle augmente les impacts
112	c radiatifs des nuages.
113	ccc LOGICAL new_oliq
114	ccc PARAMETER (new_oliq=.FALSE.)
115	c======================================================================
116	c Clefs controlant deux parametrisations de l'orographie:
117	cc LOGICAL ok_orodr
118	ccc PARAMETER (ok_orodr=.FALSE.)
119	ccc LOGICAL ok_orolf
120	ccc PARAMETER (ok_orolf=.FALSE.)
121	c======================================================================
122	LOGICAL ok_journe ! sortir le fichier journalier
123	save ok_journe
124	c PARAMETER (ok_journe=.true.)
125
126	integer lev_histday
127	save lev_histday
128	data lev_histday/1/
129	c
130	LOGICAL ok_mensuel ! sortir le fichier mensuel
131	save ok_mensuel
132	c PARAMETER (ok_mensuel=.true.)
133	c
134	LOGICAL ok_instan ! sortir le fichier instantane
135	save ok_instan
136	c PARAMETER (ok_instan=.true.)
137	c
138	LOGICAL ok_region ! sortir le fichier regional
139	PARAMETER (ok_region=.FALSE.)
140	c======================================================================
141	c
142	INTEGER ivap ! indice de traceurs pour vapeur d'eau
143	PARAMETER (ivap=1)
144	INTEGER iliq ! indice de traceurs pour eau liquide
145	PARAMETER (iliq=2)
146
147	c
148	c
149	c Variables argument:
150	c
151	INTEGER nlon
152	INTEGER nlev
153	INTEGER nqmax
154	REAL rjourvrai, rjour_ecri
155	REAL gmtime
156	REAL pdtphys
157	LOGICAL debut, lafin
158	REAL paprs(klon,klev+1)
159	REAL pplay(klon,klev)
160	REAL pphi(klon,klev)
161	REAL pphis(klon)
162	REAL paire(klon)
163	REAL presnivs(klev)
164	REAL znivsig(klev)
165	REAL zsurf(nbsrf)
166	real cufi(klon), cvfi(klon)
167
168	REAL u(klon,klev)
169	REAL v(klon,klev)
170	REAL t(klon,klev)
171	REAL qx(klon,klev,nqmax)
172
173	REAL t_ancien(klon,klev), q_ancien(klon,klev)
174	SAVE t_ancien, q_ancien
175	LOGICAL ancien_ok
176	SAVE ancien_ok
177
178	REAL d_t_dyn(klon,klev)
179	REAL d_q_dyn(klon,klev)
180
181	REAL omega(klon,klev)
182
183	REAL d_u(klon,klev)
184	REAL d_v(klon,klev)
185	REAL d_t(klon,klev)
186	REAL d_qx(klon,klev,nqmax)
187	REAL d_ps(klon)
188
189	cccIM
190	INTEGER klevp1
191	PARAMETER(klevp1=klev+1)
192	#include "raddim.h"
193	REAL*8 ZFSUP(KDLON,KFLEV+1)
194	REAL*8 ZFSDN(KDLON,KFLEV+1)
195	REAL*8 ZFSUP0(KDLON,KFLEV+1)
196	REAL*8 ZFSDN0(KDLON,KFLEV+1)
197
198	cccIM cf. FH
199	real u850(klon),v850(klon),u200(klon),v200(klon)
200	real u500(klon),v500(klon),phi500(klon),w500(klon)
201
202	logical ok_hf
203	real ecrit_hf
204	integer nid_hf
205	save ok_hf, ecrit_hf, nid_hf
206
207	c QUESTION : noms de variables ?
208
209	#define histhf
210	#ifdef histhf
211	data ok_hf,ecrit_hf/.true.,0.25/
212	#else
213	data ok_hf/.false./
214	#endif
215
216	INTEGER longcles
217	PARAMETER ( longcles = 20 )
218	REAL clesphy0( longcles )
219	c
220	c Variables quasi-arguments
221	c
222	REAL xjour
223	SAVE xjour
224	c
225	c
226	c Variables propres a la physique
227	c
228	REAL dtime
229	SAVE dtime ! pas temporel de la physique
230	c
231	INTEGER radpas
232	SAVE radpas ! frequence d'appel rayonnement
233	c
234	REAL radsol(klon)
235	SAVE radsol ! bilan radiatif au sol
236	c
237	REAL rlat(klon)
238	SAVE rlat ! latitude pour chaque point
239	c
240	REAL rlon(klon)
241	SAVE rlon ! longitude pour chaque point
242	c
243	cc INTEGER iflag_con
244	cc SAVE iflag_con ! indicateur de la convection
245	c
246	INTEGER itap
247	SAVE itap ! compteur pour la physique
248	c
249	REAL co2_ppm_etat0
250	c
251	REAL solaire_etat0
252	c
253	real slp(klon) ! sea level pressure
254
255	REAL ftsol(klon,nbsrf)
256	SAVE ftsol ! temperature du sol
257	c
258	REAL ftsoil(klon,nsoilmx,nbsrf)
259	SAVE ftsoil ! temperature dans le sol
260	c
261	REAL fevap(klon,nbsrf)
262	SAVE fevap ! evaporation
263	REAL fluxlat(klon,nbsrf)
264	SAVE fluxlat
265	c
266	REAL deltat(klon)
267	SAVE deltat ! ecart avec la SST de reference
268	c
269	REAL fqsol(klon,nbsrf)
270	SAVE fqsol ! humidite du sol
271	c
272	REAL fsnow(klon,nbsrf)
273	SAVE fsnow ! epaisseur neigeuse
274	c
275	REAL falbe(klon,nbsrf)
276	SAVE falbe ! albedo par type de surface
277	REAL falblw(klon,nbsrf)
278	SAVE falblw ! albedo par type de surface
279
280	c
281	c
282	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
283	c
284	REAL zmea(klon)
285	SAVE zmea ! orographie moyenne
286	c
287	REAL zstd(klon)
288	SAVE zstd ! deviation standard de l'OESM
289	c
290	REAL zsig(klon)
291	SAVE zsig ! pente de l'OESM
292	c
293	REAL zgam(klon)
294	save zgam ! anisotropie de l'OESM
295	c
296	REAL zthe(klon)
297	SAVE zthe ! orientation de l'OESM
298	c
299	REAL zpic(klon)
300	SAVE zpic ! Maximum de l'OESM
301	c
302	REAL zval(klon)
303	SAVE zval ! Minimum de l'OESM
304	c
305	REAL rugoro(klon)
306	SAVE rugoro ! longueur de rugosite de l'OESM
307	c
308	REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
309	c
310	REAL zuthe(klon),zvthe(klon)
311	SAVE zuthe
312	SAVE zvthe
313	INTEGER igwd,idx(klon),itest(klon)
314	c
315	REAL agesno(klon,nbsrf)
316	SAVE agesno ! age de la neige
317	c
318	REAL alb_neig(klon)
319	SAVE alb_neig ! albedo de la neige
320	cKE43
321	c Variables liees a la convection de K. Emanuel (sb):
322	c
323	REAL ema_workcbmf(klon) ! cloud base mass flux
324	SAVE ema_workcbmf
325
326	REAL ema_cbmf(klon) ! cloud base mass flux
327	SAVE ema_cbmf
328
329	REAL ema_pcb(klon) ! cloud base pressure
330	SAVE ema_pcb
331
332	REAL ema_pct(klon) ! cloud top pressure
333	SAVE ema_pct
334
335	REAL bas, top ! cloud base and top levels
336	SAVE bas
337	SAVE top
338
339	REAL Ma(klon,klev) ! undilute upward mass flux
340	SAVE Ma
341	REAL qcondc(klon,klev) ! in-cld water content from convect
342	SAVE qcondc
343	REAL ema_work1(klon, klev), ema_work2(klon, klev)
344	SAVE ema_work1, ema_work2
345	REAL wdn(klon), tdn(klon), qdn(klon)
346
347	REAL wd(klon) ! sb
348	SAVE wd ! sb
349
350	c Variables locales pour la couche limite (al1):
351	c
352	cAl1 REAL pblh(klon) ! Hauteur de couche limite
353	cAl1 SAVE pblh
354	c34EK
355	c
356	c Variables locales:
357	c
358	REAL cdragh(klon) ! drag coefficient pour T and Q
359	REAL cdragm(klon) ! drag coefficient pour vent
360	cAA
361	cAA Pour phytrac
362	cAA
363	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
364	REAL yu1(klon) ! vents dans la premiere couche U
365	REAL yv1(klon) ! vents dans la premiere couche V
366	LOGICAL offline ! Controle du stockage ds "physique"
367	PARAMETER (offline=.false.)
368	INTEGER physid
369	REAL pfrac_impa(klon,klev)! Produits des coefs lessivage impaction
370	save pfrac_impa
371	REAL pfrac_nucl(klon,klev)! Produits des coefs lessivage nucleation
372	save pfrac_nucl
373	REAL pfrac_1nucl(klon,klev)! Produits des coefs lessi nucl (alpha = 1)
374	save pfrac_1nucl
375	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
376	REAL frac_nucl(klon,klev) ! idem (nucleation)
377	cAA
378	REAL rain_fall(klon) ! pluie
379	REAL snow_fall(klon) ! neige
380	save snow_fall, rain_fall
381	REAL evap(klon), devap(klon) ! evaporation et sa derivee
382	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
383	REAL dlw(klon) ! derivee infra rouge
384	REAL bils(klon) ! bilan de chaleur au sol
385	cIM cf. JLD
386	REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
387	C type de sous-surface et pondere par la fraction
388	REAL fder(klon) ! Derive de flux (sensible et latente)
389	save fder
390	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
391	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
392	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
393	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
394	c
395	REAL frugs(klon,nbsrf) ! longueur de rugosite
396	save frugs
397	REAL zxrugs(klon) ! longueur de rugosite
398	c
399	c Conditions aux limites
400	c
401	INTEGER julien
402	c
403	INTEGER lmt_pas
404	SAVE lmt_pas ! frequence de mise a jour
405	REAL pctsrf(klon,nbsrf)
406	SAVE pctsrf ! sous-fraction du sol
407	REAL albsol(klon)
408	SAVE albsol ! albedo du sol total
409	REAL albsollw(klon)
410	SAVE albsollw ! albedo du sol total
411	REAL albsol1(klon)
412	SAVE albsol1 ! albedo du sol total
413	REAL albsollw1(klon)
414	SAVE albsollw1 ! albedo du sol total
415
416	REAL wo(klon,klev)
417	SAVE wo ! ozone
418	c======================================================================
419	c
420	c Declaration des procedures appelees
421	c
422	EXTERNAL angle ! calculer angle zenithal du soleil
423	EXTERNAL alboc ! calculer l'albedo sur ocean
424	EXTERNAL albsno ! calculer albedo sur neige
425	EXTERNAL ajsec ! ajustement sec
426	EXTERNAL clmain ! couche limite
427	EXTERNAL condsurf ! lire les conditions aux limites
428	EXTERNAL conlmd ! convection (schema LMD)
429	cKE43
430	EXTERNAL conema3 ! convect4.3
431	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
432	cAA
433	EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
434	c ! stockage des coefficients necessaires au
435	c ! lessivage OFF-LINE et ON-LINE
436	cAA
437	EXTERNAL hgardfou ! verifier les temperatures
438	EXTERNAL nuage ! calculer les proprietes radiatives
439	EXTERNAL o3cm ! initialiser l'ozone
440	EXTERNAL orbite ! calculer l'orbite terrestre
441	EXTERNAL ozonecm ! prescrire l'ozone
442	EXTERNAL phyetat0 ! lire l'etat initial de la physique
443	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
444	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
445	EXTERNAL suphec ! initialiser certaines constantes
446	EXTERNAL transp ! transport total de l'eau et de l'energie
447	EXTERNAL ecribina ! ecrire le fichier binaire global
448	EXTERNAL ecribins ! ecrire le fichier binaire global
449	EXTERNAL ecrirega ! ecrire le fichier binaire regional
450	EXTERNAL ecriregs ! ecrire le fichier binaire regional
451	c
452	c Variables locales
453	c
454	real clwcon(klon,klev),rnebcon(klon,klev)
455	real clwcon0(klon,klev),rnebcon0(klon,klev)
456	save rnebcon, clwcon
457
458	REAL rhcl(klon,klev) ! humiditi relative ciel clair
459	REAL dialiq(klon,klev) ! eau liquide nuageuse
460	REAL diafra(klon,klev) ! fraction nuageuse
461	REAL cldliq(klon,klev) ! eau liquide nuageuse
462	REAL cldfra(klon,klev) ! fraction nuageuse
463	REAL cldtau(klon,klev) ! epaisseur optique
464	REAL cldemi(klon,klev) ! emissivite infrarouge
465	c
466	CXXX PB
467	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
468	REAL fluxt(klon,klev, nbsrf) ! flux turbulent de chaleur
469	REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u
470	REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v
471	c
472	REAL zxfluxt(klon, klev)
473	REAL zxfluxq(klon, klev)
474	REAL zxfluxu(klon, klev)
475	REAL zxfluxv(klon, klev)
476	CXXX
477	REAL heat(klon,klev) ! chauffage solaire
478	REAL heat0(klon,klev) ! chauffage solaire ciel clair
479	REAL cool(klon,klev) ! refroidissement infrarouge
480	REAL cool0(klon,klev) ! refroidissement infrarouge ciel clair
481	REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)
482	real sollwdown(klon) ! downward LW flux at surface
483	REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
484	REAL albpla(klon)
485	cIM cf. JLD
486	REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
487	REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
488	c Le rayonnement n'est pas calcule tous les pas, il faut donc
489	c sauvegarder les sorties du rayonnement
490	SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
491	SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
492	cccIM
493	SAVE ZFSUP,ZFSDN,ZFSUP0,ZFSDN0
494
495	INTEGER itaprad
496	SAVE itaprad
497	c
498	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
499	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
500	c
501	REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
502	REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
503	c
504	REAL zxtsol(klon), zxqsol(klon), zxsnow(klon), zxfluxlat(klon)
505	c
506	REAL dist, rmu0(klon), fract(klon)
507	REAL zdtime, zlongi
508	c
509	CHARACTER*2 str2
510	CHARACTER*2 iqn
511	c
512	REAL qcheck
513	REAL z_avant(klon), z_apres(klon), z_factor(klon)
514	LOGICAL zx_ajustq
515	c
516	REAL za, zb
517	REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
518	real zqsat(klon,klev)
519	INTEGER i, k, iq, ig, j, nsrf, ll
520	REAL t_coup
521	PARAMETER (t_coup=234.0)
522	c
523	REAL zphi(klon,klev)
524	REAL zx_tmp_x(iim), zx_tmp_yjjmp1
525	REAL zx_relief(iim,jjmp1)
526	REAL zx_aire(iim,jjmp1)
527	cKE43
528	c Variables locales pour la convection de K. Emanuel (sb):
529	c
530	REAL upwd(klon,klev) ! saturated updraft mass flux
531	REAL dnwd(klon,klev) ! saturated downdraft mass flux
532	REAL dnwd0(klon,klev) ! unsaturated downdraft mass flux
533	REAL tvp(klon,klev) ! virtual temp of lifted parcel
534	REAL cape(klon) ! CAPE
535	SAVE cape
536	cccIM
537	CHARACTER*40 capemaxcels
538
539	REAL pbase(klon) ! cloud base pressure
540	SAVE pbase
541	REAL bbase(klon) ! cloud base buoyancy
542	SAVE bbase
543	REAL rflag(klon) ! flag fonctionnement de convect
544	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
545	c -- convect43:
546	INTEGER ntra ! nb traceurs pour convect4.3
547	REAL pori_con(klon) ! pressure at the origin level of lifted parcel
548	REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
549	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
550	REAL dplcldt(klon), dplcldr(klon)
551	c? . condm_con(klon,klev),conda_con(klon,klev),
552	c? . mr_con(klon,klev),ep_con(klon,klev)
553	c? . ,sadiab(klon,klev),wadiab(klon,klev)
554	c --
555	c34EK
556	c
557	c Variables du changement
558	c
559	c con: convection
560	c lsc: condensation a grande echelle (Large-Scale-Condensation)
561	c ajs: ajustement sec
562	c eva: evaporation de l'eau liquide nuageuse
563	c vdf: couche limite (Vertical DiFfusion)
564	REAL d_t_con(klon,klev),d_q_con(klon,klev)
565	REAL d_u_con(klon,klev),d_v_con(klon,klev)
566	REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev)
567	REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
568	REAL d_t_eva(klon,klev),d_q_eva(klon,klev)
569	REAL rneb(klon,klev)
570	c
571	REAL pmfu(klon,klev), pmfd(klon,klev)
572	REAL pen_u(klon,klev), pen_d(klon,klev)
573	REAL pde_u(klon,klev), pde_d(klon,klev)
574	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
575	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
576	REAL prfl(klon,klev+1), psfl(klon,klev+1)
577	c
578	INTEGER ibas_con(klon), itop_con(klon)
579	REAL rain_con(klon), rain_lsc(klon)
580	REAL snow_con(klon), snow_lsc(klon)
581	REAL d_ts(klon,nbsrf)
582	c
583	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
584	REAL d_t_vdf(klon,klev), d_q_vdf(klon,klev)
585	c
586	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
587	REAL d_t_oro(klon,klev)
588	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
589	REAL d_t_lif(klon,klev)
590
591	REAL ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev)
592	real ratqsbas,ratqshaut
593	save ratqsbas,ratqshaut, ratqs
594	real zpt_conv(klon,klev)
595
596	c Parametres lies au nouveau schema de nuages (SB, PDF)
597	real fact_cldcon
598	real facttemps
599	logical ok_newmicro
600	save ok_newmicro
601	save fact_cldcon,facttemps
602	real facteur
603
604	integer iflag_cldcon
605	save iflag_cldcon
606
607	logical ptconv(klon,klev)
608
609	c
610	c Variables liees a l'ecriture de la bande histoire physique
611	c
612	INTEGER ecrit_mth
613	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
614	c
615	INTEGER ecrit_day
616	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
617	c
618	INTEGER ecrit_ins
619	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
620	c
621	INTEGER ecrit_reg
622	SAVE ecrit_reg ! frequence d'ecriture
623	c
624	integer itau_w ! pas de temps ecriture = itap + itau_phy
625	c
626	c
627	c Variables locales pour effectuer les appels en serie
628	c
629	REAL t_seri(klon,klev), q_seri(klon,klev)
630	REAL ql_seri(klon,klev),qs_seri(klon,klev)
631	REAL u_seri(klon,klev), v_seri(klon,klev)
632	c
633	REAL tr_seri(klon,klev,nbtr)
634	REAL d_tr(klon,klev,nbtr)
635
636	REAL zx_rh(klon,klev)
637
638	INTEGER length
639	PARAMETER ( length = 100 )
640	REAL tabcntr0( length )
641	c
642	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
643	REAL zx_tmp_fi2d(klon)
644	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
645	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
646	c
647	INTEGER nid_day, nid_mth, nid_ins
648	SAVE nid_day, nid_mth, nid_ins
649	c
650	INTEGER nhori, nvert
651	REAL zsto, zout
652	real zjulian
653	save zjulian
654
655	character*20 modname
656	character*80 abort_message
657	logical ok_sync
658	real date0
659	integer idayref
660
661	C essai writephys
662	integer fid_day, fid_mth, fid_ins
663	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
664	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
665	parameter (prof2d_on = 1, prof3d_on = 2,
666	. prof2d_av = 3, prof3d_av = 4)
667	character*30 nom_fichier
668	character*10 varname
669	character*40 vartitle
670	character*20 varunits
671	C Variables liees au bilan d'energie et d'enthalpi
672	REAL ztsol(klon)
673	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
674	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
675	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
676	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
677	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
678	REAL d_h_vcol_phy
679	REAL fs_bound, fq_bound
680	SAVE d_h_vcol_phy
681	REAL zero_v(klon)
682	CHARACTER*15 ztit
683	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
684	SAVE ip_ebil
685	DATA ip_ebil/2/
686	INTEGER if_ebil ! level for energy conserv. dignostics
687	SAVE if_ebil
688	c+jld ec_conser
689	REAL d_t_ec(klon,klev) ! tendance du a la conersion Ec -> E thermique
690	REAL ZRCPD
691	c-jld ec_conser
692	cIM
693	REAL t2m(klon,nbsrf), q2m(klon,nbsrf)
694	REAL u10m(klon,nbsrf), v10m(klon,nbsrf)
695	REAL zt2m(klon), zq2m(klon)
696	REAL zu10m(klon), zv10m(klon)
697	CHARACTER*40 t2mincels, t2maxcels
698	c
699	c Declaration des constantes et des fonctions thermodynamiques
700	c
701	#include "YOMCST.h"
702	#include "YOETHF.h"
703	#include "FCTTRE.h"
704	c======================================================================
705	modname = 'physiq'
706	IF (if_ebil.ge.1) THEN
707	DO i=1,klon
708	zero_v(i)=0.
709	END DO
710	END IF
711	ok_sync=.TRUE.
712	IF (nqmax .LT. 2) THEN
713	PRINT*, 'eaux vapeur et liquide sont indispensables'
714	CALL ABORT
715	ENDIF
716	IF (debut) THEN
717	CALL suphec ! initialiser constantes et parametres phys.
718	ENDIF
719
720
721	c======================================================================
722	xjour = rjourvrai
723	c
724	c Si c'est le debut, il faut initialiser plusieurs choses
725	c ********
726	c
727	IF (debut) THEN
728	C
729	IF (if_ebil.ge.1) d_h_vcol_phy=0.
730	c
731	c appel a la lecture du run.def physique
732	c
733	call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
734	. ok_instan, fact_cldcon, facttemps,ok_newmicro,
735	. iflag_cldcon,ratqsbas,ratqshaut, if_ebil)
736	cIM . , RI0)
737
738	c
739	c
740	c Initialiser les compteurs:
741	c
742
743	frugs = 0.
744	itap = 0
745	itaprad = 0
746	c
747	CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
748	. rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsol,fsnow,
749	. falbe, fevap, rain_fall,snow_fall,solsw, sollwdown,
750	. dlw,radsol,frugs,agesno,clesphy0,
751	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
752	. t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon )
753
754	c
755	radpas = NINT( 86400./dtime/nbapp_rad)
756
757	c
758	CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe,
759	, ok_instan, ok_region )
760	c
761	IF (ABS(dtime-pdtphys).GT.0.001) THEN
762	PRINT*, 'Pas physique n est pas correcte',dtime,pdtphys
763	abort_message=' See above '
764	call abort_gcm(modname,abort_message,1)
765	ENDIF
766	IF (nlon .NE. klon) THEN
767	PRINT*, 'nlon et klon ne sont pas coherents', nlon, klon
768	abort_message=' See above '
769	call abort_gcm(modname,abort_message,1)
770	ENDIF
771	IF (nlev .NE. klev) THEN
772	PRINT*, 'nlev et klev ne sont pas coherents', nlev, klev
773	abort_message=' See above '
774	call abort_gcm(modname,abort_message,1)
775	ENDIF
776	c
777	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
778	PRINT*, 'Nbre d appels au rayonnement insuffisant'
779	PRINT*, "Au minimum 4 appels par jour si cycle diurne"
780	abort_message=' See above '
781	call abort_gcm(modname,abort_message,1)
782	ENDIF
783	PRINT*, "Clef pour la convection, iflag_con=", iflag_con
784	PRINT*, "Clef pour le driver de la convection, ok_cvl=", ok_cvl
785	c
786	cKE43
787	c Initialisation pour la convection de K.E. (sb):
788	IF (iflag_con.GE.3) THEN
789
790	PRINT, "** Convection de Kerry Emanuel 4.3 "
791	PRINT*, "On va utiliser le melange convectif des traceurs qui"
792	PRINT*, "est calcule dans convect4.3"
793	PRINT*, " !!! penser aux logical flags de phytrac"
794
795	DO i = 1, klon
796	ema_cbmf(i) = 0.
797	ema_pcb(i) = 0.
798	ema_pct(i) = 0.
799	ema_workcbmf(i) = 0.
800	ENDDO
801
802	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
803	DO i = 1, klon
804	ibas_con(i) = 1
805	itop_con(i) = klev+1
806	ENDDO
807	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
808
809	ENDIF
810
811	c34EK
812	IF (ok_orodr) THEN
813	DO i=1,klon
814	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
815	ENDDO
816	CALL SUGWD(klon,klev,paprs,pplay)
817	DO i=1,klon
818	zuthe(i)=0.
819	zvthe(i)=0.
820	if(zstd(i).gt.10.)then
821	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
822	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
823	endif
824	ENDDO
825	ENDIF
826	c
827	c
828	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
829	PRINT*,'La frequence de lecture surface est de ', lmt_pas
830	c
831	ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
832	IF (ok_mensuel) THEN
833	PRINT*, 'La frequence de sortie mensuelle est de ', ecrit_mth
834	ENDIF
835	ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
836	IF (ok_journe) THEN
837	PRINT*, 'La frequence de sortie journaliere est de ',ecrit_day
838	ENDIF
839	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
840	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
841	ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
842	ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
843	IF (ok_instan) THEN
844	PRINT*, 'La frequence de sortie instant. est de ', ecrit_ins
845	ENDIF
846	ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
847	IF (ok_region) THEN
848	PRINT*, 'La frequence de sortie region est de ', ecrit_reg
849	ENDIF
850
851	c
852	c Initialiser le couplage si necessaire
853	c
854	npas = 0
855	nexca = 0
856	if (ocean == 'couple') then
857	npas = itaufin/ iphysiq
858	nexca = 86400 / dtime
859	write(,)' ##### Ocean couple #####'
860	write(,)' Valeurs des pas de temps'
861	write(,)' npas = ', npas
862	write(,)' nexca = ', nexca
863	endif
864	c
865	c
866	cccIM
867	capemaxcels = 't_max(X)'
868	t2mincels = 't_min(X)'
869	t2maxcels = 't_max(X)'
870
871	cccIM cf. FH
872	c
873	c=============================================================
874	c Initialisation des sorties
875	c=============================================================
876	#ifdef histhf
877	#include "ini_histhf.h"
878	#endif
879
880	#include "ini_histday.h"
881	#include "ini_histmth.h"
882	#include "ini_histins.h"
883
884	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
885	date0 = zjulian
886	C date0 = day_ini
887	WRITE(,) 'physiq date0 : ',date0
888	c
889	c
890	c
891	c Prescrire l'ozone dans l'atmosphere
892	c
893	c
894	cc DO i = 1, klon
895	cc DO k = 1, klev
896	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
897	cc ENDDO
898	cc ENDDO
899	c
900	c
901	ENDIF
902	c
903	c ************** Fin de IF ( debut ) *************
904	c
905	c
906	c Mettre a zero des variables de sortie (pour securite)
907	c
908	DO i = 1, klon
909	d_ps(i) = 0.0
910	ENDDO
911	DO k = 1, klev
912	DO i = 1, klon
913	d_t(i,k) = 0.0
914	d_u(i,k) = 0.0
915	d_v(i,k) = 0.0
916	ENDDO
917	ENDDO
918	DO iq = 1, nqmax
919	DO k = 1, klev
920	DO i = 1, klon
921	d_qx(i,k,iq) = 0.0
922	ENDDO
923	ENDDO
924	ENDDO
925	c
926	c Ne pas affecter les valeurs entrees de u, v, h, et q
927	c
928	DO k = 1, klev
929	DO i = 1, klon
930	t_seri(i,k) = t(i,k)
931	u_seri(i,k) = u(i,k)
932	v_seri(i,k) = v(i,k)
933	q_seri(i,k) = qx(i,k,ivap)
934	ql_seri(i,k) = qx(i,k,iliq)
935	qs_seri(i,k) = 0.
936	ENDDO
937	ENDDO
938	IF (nqmax.GE.3) THEN
939	DO iq = 3, nqmax
940	DO k = 1, klev
941	DO i = 1, klon
942	tr_seri(i,k,iq-2) = qx(i,k,iq)
943	ENDDO
944	ENDDO
945	ENDDO
946	ELSE
947	DO k = 1, klev
948	DO i = 1, klon
949	tr_seri(i,k,1) = 0.0
950	ENDDO
951	ENDDO
952	ENDIF
953	C
954	IF (if_ebil.ge.1) THEN
955	DO i = 1, klon
956	ztsol(i) = 0.
957	ENDDO
958	DO nsrf = 1, nbsrf
959	DO i = 1, klon
960	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
961	ENDDO
962	ENDDO
963	ztit='after dynamic'
964	CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
965	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
966	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
967	C Comme les tendances de la physique sont ajoute dans la dynamique,
968	C on devrait avoir que la variation d'entalpie par la dynamique
969	C est egale a la variation de la physique au pas de temps precedent.
970	C Donc la somme de ces 2 variations devrait etre nulle.
971	call diagphy(paire,ztit,ip_ebil
972	e , zero_v, zero_v, zero_v, zero_v, zero_v
973	e , zero_v, zero_v, zero_v, ztsol
974	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
975	s , fs_bound, fq_bound )
976	END IF
977
978	c Diagnostiquer la tendance dynamique
979	c
980	IF (ancien_ok) THEN
981	DO k = 1, klev
982	DO i = 1, klon
983	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
984	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
985	ENDDO
986	ENDDO
987	ELSE
988	DO k = 1, klev
989	DO i = 1, klon
990	d_t_dyn(i,k) = 0.0
991	d_q_dyn(i,k) = 0.0
992	ENDDO
993	ENDDO
994	ancien_ok = .TRUE.
995	ENDIF
996	c
997	c Ajouter le geopotentiel du sol:
998	c
999	DO k = 1, klev
1000	DO i = 1, klon
1001	zphi(i,k) = pphi(i,k) + pphis(i)
1002	ENDDO
1003	ENDDO
1004	c
1005	c Verifier les temperatures
1006	c
1007	CALL hgardfou(t_seri,ftsol,'debutphy')
1008	c
1009	c Incrementer le compteur de la physique
1010	c
1011	itap = itap + 1
1012	julien = MOD(NINT(xjour),360)
1013	c
1014	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1015	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1016	c
1017	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1018	PRINT *,' PHYS cond julien ',julien
1019	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1020	ENDIF
1021	c
1022	c Re-evaporer l'eau liquide nuageuse
1023	c
1024	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1025	DO i = 1, klon
1026	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1027	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1028	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1029	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1030	zb = MAX(0.0,ql_seri(i,k))
1031	za = - MAX(0.0,ql_seri(i,k))
1032	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1033	t_seri(i,k) = t_seri(i,k) + za
1034	q_seri(i,k) = q_seri(i,k) + zb
1035	ql_seri(i,k) = 0.0
1036	d_t_eva(i,k) = za
1037	d_q_eva(i,k) = zb
1038	ENDDO
1039	ENDDO
1040	c
1041	IF (if_ebil.ge.2) THEN
1042	ztit='after reevap'
1043	CALL diagetpq(paire,ztit,ip_ebil,2,1,dtime
1044	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1045	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1046	call diagphy(paire,ztit,ip_ebil
1047	e , zero_v, zero_v, zero_v, zero_v, zero_v
1048	e , zero_v, zero_v, zero_v, ztsol
1049	e , d_h_vcol, d_qt, d_ec
1050	s , fs_bound, fq_bound )
1051	C
1052	END IF
1053	C
1054	c
1055	c Appeler la diffusion verticale (programme de couche limite)
1056	c
1057	DO i = 1, klon
1058	c if (.not. ok_veget) then
1059	c frugs(i,is_ter) = SQRT(frugs(i,is_ter)2+rugoro(i)2)
1060	c endif
1061	c frugs(i,is_lic) = rugoro(i)
1062	c frugs(i,is_oce) = rugmer(i)
1063	c frugs(i,is_sic) = 0.001
1064	zxrugs(i) = 0.0
1065	ENDDO
1066	DO nsrf = 1, nbsrf
1067	DO i = 1, klon
1068	frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
1069	cccc frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
1070	ENDDO
1071	ENDDO
1072	DO nsrf = 1, nbsrf
1073	DO i = 1, klon
1074	zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
1075	ENDDO
1076	ENDDO
1077	c
1078	C calculs necessaires au calcul de l'albedo dans l'interface
1079	c
1080	CALL orbite(FLOAT(julien),zlongi,dist)
1081	IF (cycle_diurne) THEN
1082	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1083	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1084	ELSE
1085	rmu0 = -999.999
1086	ENDIF
1087
1088	DO nsrf = 1, nbsrf
1089	DO i = 1, klon
1090	fsollw(i,nsrf) = sollwdown(i) - RSIGMAftsol(i,nsrf)*4
1091	fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
1092	ENDDO
1093	ENDDO
1094
1095	fder = dlw
1096
1097	CALL clmain(dtime,itap,date0,pctsrf,
1098	e t_seri,q_seri,u_seri,v_seri,
1099	e julien, rmu0,
1100	e ok_veget, ocean, npas, nexca, ftsol,
1101	$ soil_model,ftsoil,
1102	$ paprs,pplay,radsol, fsnow,fqsol,fevap,falbe,falblw,
1103	$ fluxlat,
1104	cIM cf. JLD e rain_fall, snow_fall, solsw, sollw, sollwdown, fder,
1105	e rain_fall, snow_fall, fsolsw, fsollw, sollwdown, fder,
1106	e rlon, rlat, cufi, cvfi, frugs,
1107	e debut, lafin, agesno,rugoro ,
1108	s d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
1109	s fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
1110	s dsens, devap,
1111	s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m)
1112	c
1113	CXXX PB
1114	CXXX Incrementation des flux
1115	CXXX
1116	zxfluxt=0.
1117	zxfluxq=0.
1118	zxfluxu=0.
1119	zxfluxv=0.
1120	DO nsrf = 1, nbsrf
1121	DO k = 1, klev
1122	DO i = 1, klon
1123	zxfluxt(i,k) = zxfluxt(i,k) +
1124	$ fluxt(i,k,nsrf) * pctsrf( i, nsrf)
1125	zxfluxq(i,k) = zxfluxq(i,k) +
1126	$ fluxq(i,k,nsrf) * pctsrf( i, nsrf)
1127	zxfluxu(i,k) = zxfluxu(i,k) +
1128	$ fluxu(i,k,nsrf) * pctsrf( i, nsrf)
1129	zxfluxv(i,k) = zxfluxv(i,k) +
1130	$ fluxv(i,k,nsrf) * pctsrf( i, nsrf)
1131	END DO
1132	END DO
1133	END DO
1134	DO i = 1, klon
1135	sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1136	c evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
1137	evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1138	fder(i) = dlw(i) + dsens(i) + devap(i)
1139	ENDDO
1140
1141
1142	DO k = 1, klev
1143	DO i = 1, klon
1144	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1145	q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
1146	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1147	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1148	ENDDO
1149	ENDDO
1150	c
1151	IF (if_ebil.ge.2) THEN
1152	ztit='after clmain'
1153	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1154	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1155	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1156	call diagphy(paire,ztit,ip_ebil
1157	e , zero_v, zero_v, zero_v, zero_v, sens
1158	e , evap , zero_v, zero_v, ztsol
1159	e , d_h_vcol, d_qt, d_ec
1160	s , fs_bound, fq_bound )
1161	END IF
1162	C
1163	c
1164	c Incrementer la temperature du sol
1165	c
1166	DO i = 1, klon
1167	zxtsol(i) = 0.0
1168	zxfluxlat(i) = 0.0
1169	cccIM
1170	zt2m(i) = 0.0
1171	zq2m(i) = 0.0
1172	zu10m(i) = 0.0
1173	zv10m(i) = 0.0
1174	c
1175	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
1176	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
1177	$ THEN
1178	WRITE(,) 'physiq : pb sous surface au point ', i,
1179	$ pctsrf(i, 1 : nbsrf)
1180	ENDIF
1181	ENDDO
1182	DO nsrf = 1, nbsrf
1183	DO i = 1, klon
1184	c IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1185	ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
1186	cIM cf. JLD
1187	wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf)
1188	$ + fluxt(i,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1189	zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1190	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf)
1191	cccIM
1192	zt2m(i) = zt2m(i) + t2m(i,nsrf)*pctsrf(i,nsrf)
1193	zq2m(i) = zq2m(i) + q2m(i,nsrf)*pctsrf(i,nsrf)
1194	zu10m(i) = zu10m(i) + u10m(i,nsrf)*pctsrf(i,nsrf)
1195	zv10m(i) = zv10m(i) + v10m(i,nsrf)*pctsrf(i,nsrf)
1196	c ENDIF
1197	ENDDO
1198	ENDDO
1199
1200	c
1201	c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
1202	c
1203	DO nsrf = 1, nbsrf
1204	DO i = 1, klon
1205	IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i)
1206	cccIM
1207	IF (pctsrf(i,nsrf) .LT. epsfra) t2m(i,nsrf) = zt2m(i)
1208	IF (pctsrf(i,nsrf) .LT. epsfra) q2m(i,nsrf) = zq2m(i)
1209	IF (pctsrf(i,nsrf) .LT. epsfra) u10m(i,nsrf) = zu10m(i)
1210	IF (pctsrf(i,nsrf) .LT. epsfra) v10m(i,nsrf) = zv10m(i)
1211	ENDDO
1212	ENDDO
1213	c
1214	c
1215	c Calculer la derive du flux infrarouge
1216	c
1217	cXXX DO nsrf = 1, nbsrf
1218	DO i = 1, klon
1219	cXXX IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1220	dlw(i) = - 4.0RSIGMAzxtsol(i)**3
1221	cXXX . *(ftsol(i,nsrf)-zxtsol(i))
1222	cXXX . *pctsrf(i,nsrf)
1223	cXXX ENDIF
1224	cXXX ENDDO
1225	ENDDO
1226	c
1227	c Appeler la convection (au choix)
1228	c
1229	DO k = 1, klev
1230	DO i = 1, klon
1231	conv_q(i,k) = d_q_dyn(i,k)
1232	. + d_q_vdf(i,k)/dtime
1233	conv_t(i,k) = d_t_dyn(i,k)
1234	. + d_t_vdf(i,k)/dtime
1235	ENDDO
1236	ENDDO
1237	IF (check) THEN
1238	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1239	PRINT*, "avantcon=", za
1240	ENDIF
1241	zx_ajustq = .FALSE.
1242	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1243	IF (zx_ajustq) THEN
1244	DO i = 1, klon
1245	z_avant(i) = 0.0
1246	ENDDO
1247	DO k = 1, klev
1248	DO i = 1, klon
1249	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1250	. *(paprs(i,k)-paprs(i,k+1))/RG
1251	ENDDO
1252	ENDDO
1253	ENDIF
1254	IF (iflag_con.EQ.1) THEN
1255	stop'reactiver le call conlmd dans physiq.F'
1256	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1257	c . d_t_con, d_q_con,
1258	c . rain_con, snow_con, ibas_con, itop_con)
1259	ELSE IF (iflag_con.EQ.2) THEN
1260	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1261	e conv_t, conv_q, zxfluxq(1,1), omega,
1262	s d_t_con, d_q_con, rain_con, snow_con,
1263	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1264	s kcbot, kctop, kdtop, pmflxr, pmflxs)
1265	WHERE (rain_con < 0.) rain_con = 0.
1266	WHERE (snow_con < 0.) snow_con = 0.
1267	DO i = 1, klon
1268	ibas_con(i) = klev+1 - kcbot(i)
1269	itop_con(i) = klev+1 - kctop(i)
1270	ENDDO
1271	ELSE IF (iflag_con.GE.3) THEN
1272	c nb of tracers for the KE convection:
1273	if (nqmax .GE. 4) then
1274	ntra = nbtr
1275	else
1276	ntra = 1
1277	endif
1278	c
1279	c sb, oct02:
1280	c Schema de convection modularise et vectorise:
1281	c (driver commun aux versions 3 et 4)
1282	c
1283	IF (ok_cvl) THEN ! new driver for convectL
1284
1285	CALL concvl (iflag_con,
1286	. dtime,paprs,pplay,t_seri,q_seri,
1287	. u_seri,v_seri,tr_seri,nbtr,
1288	. ema_work1,ema_work2,
1289	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1290	. rain_con, snow_con, ibas_con, itop_con,
1291	. upwd,dnwd,dnwd0,
1292	. Ma,cape,tvp,iflagctrl,
1293	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd)
1294	cIM cf. FH
1295	clwcon0=qcondc
1296
1297	ELSE ! ok_cvl
1298
1299	c print*,'Avant conema OUI'
1300	CALL conema3 (dtime,
1301	. paprs,pplay,t_seri,q_seri,
1302	. u_seri,v_seri,tr_seri,nbtr,
1303	. ema_work1,ema_work2,
1304	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1305	. rain_con, snow_con, ibas_con, itop_con,
1306	. upwd,dnwd,dnwd0,bas,top,
1307	. Ma,cape,tvp,rflag,
1308	. pbase
1309	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
1310	. ,clwcon0)
1311	print*,'Apres conema3 '
1312
1313	ENDIF ! ok_cvl
1314
1315	IF (.NOT. ok_gust) THEN
1316	do i = 1, klon
1317	wd(i)=0.0
1318	enddo
1319	ENDIF
1320
1321	c =================================================================== c
1322	c Calcul des proprietes des nuages convectifs
1323	c
1324	DO k = 1, klev
1325	DO i = 1, klon
1326	zx_t = t_seri(i,k)
1327	IF (thermcep) THEN
1328	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1329	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1330	zx_qs = MIN(0.5,zx_qs)
1331	zcor = 1./(1.-retv*zx_qs)
1332	zx_qs = zx_qs*zcor
1333	ELSE
1334	IF (zx_t.LT.t_coup) THEN
1335	zx_qs = qsats(zx_t)/pplay(i,k)
1336	ELSE
1337	zx_qs = qsatl(zx_t)/pplay(i,k)
1338	ENDIF
1339	ENDIF
1340	zqsat(i,k)=zx_qs
1341	ENDDO
1342	ENDDO
1343
1344	c calcul des proprietes des nuages convectifs
1345	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
1346	call clouds_gno
1347	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
1348
1349	c =================================================================== c
1350
1351	DO i = 1, klon
1352	ema_pcb(i) = pbase(i)
1353	ENDDO
1354	DO i = 1, klon
1355	ema_pct(i) = paprs(i,itop_con(i))
1356	ENDDO
1357	DO i = 1, klon
1358	ema_cbmf(i) = ema_workcbmf(i)
1359	ENDDO
1360	ELSE
1361	PRINT*, "iflag_con non-prevu", iflag_con
1362	CALL abort
1363	ENDIF
1364
1365	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
1366	c . d_u_con, d_v_con)
1367	DO k = 1, klev
1368	DO i = 1, klon
1369	t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
1370	q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
1371	u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
1372	v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
1373	ENDDO
1374	ENDDO
1375	c
1376	IF (if_ebil.ge.2) THEN
1377	ztit='after convect'
1378	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1379	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1380	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1381	call diagphy(paire,ztit,ip_ebil
1382	e , zero_v, zero_v, zero_v, zero_v, zero_v
1383	e , zero_v, rain_con, snow_con, ztsol
1384	e , d_h_vcol, d_qt, d_ec
1385	s , fs_bound, fq_bound )
1386	END IF
1387	C
1388	IF (check) THEN
1389	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1390	PRINT*, "aprescon=", za
1391	zx_t = 0.0
1392	za = 0.0
1393	DO i = 1, klon
1394	za = za + paire(i)/FLOAT(klon)
1395	zx_t = zx_t + (rain_con(i)+snow_con(i))*paire(i)/FLOAT(klon)
1396	ENDDO
1397	zx_t = zx_t/za*dtime
1398	PRINT*, "Precip=", zx_t
1399	ENDIF
1400	IF (zx_ajustq) THEN
1401	DO i = 1, klon
1402	z_apres(i) = 0.0
1403	ENDDO
1404	DO k = 1, klev
1405	DO i = 1, klon
1406	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
1407	. *(paprs(i,k)-paprs(i,k+1))/RG
1408	ENDDO
1409	ENDDO
1410	DO i = 1, klon
1411	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
1412	. /z_apres(i)
1413	ENDDO
1414	DO k = 1, klev
1415	DO i = 1, klon
1416	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
1417	. z_factor(i).LT.(1.0-1.0E-08)) THEN
1418	q_seri(i,k) = q_seri(i,k) * z_factor(i)
1419	ENDIF
1420	ENDDO
1421	ENDDO
1422	ENDIF
1423	zx_ajustq=.FALSE.
1424	c
1425	IF (nqmax.GT.2) THEN !--melange convectif de traceurs
1426	c
1427	IF (iflag_con .NE. 2 .AND. debut) THEN
1428	PRINT*, 'Pour l instant, seul conflx fonctionne ',
1429	$ 'avec traceurs', iflag_con
1430	PRINT*,' Mettre iflag_con',
1431	$ ' = 2 dans run.def et repasser'
1432	c CALL abort
1433	ENDIF
1434	c
1435	ENDIF !--nqmax.GT.2
1436	c
1437	c Appeler l'ajustement sec
1438	c
1439	CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs)
1440	DO k = 1, klev
1441	DO i = 1, klon
1442	t_seri(i,k) = t_seri(i,k) + d_t_ajs(i,k)
1443	q_seri(i,k) = q_seri(i,k) + d_q_ajs(i,k)
1444	ENDDO
1445	ENDDO
1446	c
1447	IF (if_ebil.ge.2) THEN
1448	ztit='after dry_adjust'
1449	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1450	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1451	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1452	END IF
1453
1454
1455	c-------------------------------------------------------------------------
1456	c Caclul des ratqs
1457	c-------------------------------------------------------------------------
1458
1459	c print*,'calcul des ratqs'
1460	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
1461	c ----------------
1462	c on ecrase le tableau ratqsc calcule par clouds_gno
1463	if (iflag_cldcon.eq.1) then
1464	do k=1,klev
1465	do i=1,klon
1466	if(ptconv(i,k)) then
1467	ratqsc(i,k)=ratqsbas
1468	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
1469	else
1470	ratqsc(i,k)=0.
1471	endif
1472	enddo
1473	enddo
1474	endif
1475
1476	c ratqs stables
1477	c -------------
1478	do k=1,klev
1479	ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
1480	s min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)
1481	enddo
1482
1483
1484	c ratqs final
1485	c -----------
1486	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
1487	c les ratqs sont une conbinaison de ratqss et ratqsc
1488	c ratqs final
1489	c 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
1490	c relaxation des ratqs
1491	c facttemps=exp(-pdtphys/1.e4)
1492	facteur=exp(-pdtphys*facttemps)
1493	ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
1494	ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
1495	c print*,'calcul des ratqs fini'
1496	else
1497	c on ne prend que le ratqs stable pour fisrtilp
1498	ratqs(:,:)=ratqss(:,:)
1499	endif
1500
1501
1502	c
1503	c Appeler le processus de condensation a grande echelle
1504	c et le processus de precipitation
1505	c-------------------------------------------------------------------------
1506	CALL fisrtilp(dtime,paprs,pplay,
1507	. t_seri, q_seri,ptconv,ratqs,
1508	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
1509	. rain_lsc, snow_lsc,
1510	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
1511	. frac_impa, frac_nucl,
1512	. prfl, psfl, rhcl)
1513
1514	WHERE (rain_lsc < 0) rain_lsc = 0.
1515	WHERE (snow_lsc < 0) snow_lsc = 0.
1516	DO k = 1, klev
1517	DO i = 1, klon
1518	t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
1519	q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
1520	ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
1521	cldfra(i,k) = rneb(i,k)
1522	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
1523	ENDDO
1524	ENDDO
1525	IF (check) THEN
1526	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1527	PRINT*, "apresilp=", za
1528	zx_t = 0.0
1529	za = 0.0
1530	DO i = 1, klon
1531	za = za + paire(i)/FLOAT(klon)
1532	zx_t = zx_t + (rain_lsc(i)+snow_lsc(i))*paire(i)/FLOAT(klon)
1533	ENDDO
1534	zx_t = zx_t/za*dtime
1535	PRINT*, "Precip=", zx_t
1536	ENDIF
1537	c
1538	IF (if_ebil.ge.2) THEN
1539	ztit='after fisrt'
1540	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1541	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1542	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1543	call diagphy(paire,ztit,ip_ebil
1544	e , zero_v, zero_v, zero_v, zero_v, zero_v
1545	e , zero_v, rain_lsc, snow_lsc, ztsol
1546	e , d_h_vcol, d_qt, d_ec
1547	s , fs_bound, fq_bound )
1548	END IF
1549	c
1550	c-------------------------------------------------------------------
1551	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1552	c-------------------------------------------------------------------
1553
1554	c 1. NUAGES CONVECTIFS
1555	c
1556	IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
1557
1558	c Nuages diagnostiques pour Tiedtke
1559	CALL diagcld1(paprs,pplay,
1560	. rain_con,snow_con,ibas_con,itop_con,
1561	. diafra,dialiq)
1562	DO k = 1, klev
1563	DO i = 1, klon
1564	IF (diafra(i,k).GT.cldfra(i,k)) THEN
1565	cldliq(i,k) = dialiq(i,k)
1566	cldfra(i,k) = diafra(i,k)
1567	ENDIF
1568	ENDDO
1569	ENDDO
1570
1571	ELSE IF (iflag_cldcon.eq.3) THEN
1572	c On prend pour les nuages convectifs le max du calcul de la
1573	c convection et du calcul du pas de temps précédent diminué d'un facteur
1574	c facttemps
1575	c facttemps=pdtphys/1.e4
1576	facteur = pdtphys *facttemps
1577	do k=1,klev
1578	do i=1,klon
1579	rnebcon(i,k)=rnebcon(i,k)*facteur
1580	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
1581	s then
1582	rnebcon(i,k)=rnebcon0(i,k)
1583	clwcon(i,k)=clwcon0(i,k)
1584	endif
1585	enddo
1586	enddo
1587
1588	c On prend la somme des fractions nuageuses et des contenus en eau
1589	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
1590	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
1591
1592
1593	ENDIF
1594	c
1595	c 2. NUAGES STARTIFORMES
1596	c
1597	IF (ok_stratus) THEN
1598	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
1599	DO k = 1, klev
1600	DO i = 1, klon
1601	IF (diafra(i,k).GT.cldfra(i,k)) THEN
1602	cldliq(i,k) = dialiq(i,k)
1603	cldfra(i,k) = diafra(i,k)
1604	ENDIF
1605	ENDDO
1606	ENDDO
1607	ENDIF
1608	c
1609	c Precipitation totale
1610	c
1611	DO i = 1, klon
1612	rain_fall(i) = rain_con(i) + rain_lsc(i)
1613	snow_fall(i) = snow_con(i) + snow_lsc(i)
1614	ENDDO
1615	c
1616	IF (if_ebil.ge.2) THEN
1617	ztit="after diagcld"
1618	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1619	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1620	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1621	END IF
1622	c
1623	c Calculer l'humidite relative pour diagnostique
1624	c
1625	DO k = 1, klev
1626	DO i = 1, klon
1627	zx_t = t_seri(i,k)
1628	IF (thermcep) THEN
1629	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1630	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1631	zx_qs = MIN(0.5,zx_qs)
1632	zcor = 1./(1.-retv*zx_qs)
1633	zx_qs = zx_qs*zcor
1634	ELSE
1635	IF (zx_t.LT.t_coup) THEN
1636	zx_qs = qsats(zx_t)/pplay(i,k)
1637	ELSE
1638	zx_qs = qsatl(zx_t)/pplay(i,k)
1639	ENDIF
1640	ENDIF
1641	zx_rh(i,k) = q_seri(i,k)/zx_qs
1642	zqsat(i,k)=zx_qs
1643	ENDDO
1644	ENDDO
1645	c
1646	c Calculer les parametres optiques des nuages et quelques
1647	c parametres pour diagnostiques:
1648	c
1649	if (ok_newmicro) then
1650	CALL newmicro (paprs, pplay,ok_newmicro,
1651	. t_seri, cldliq, cldfra, cldtau, cldemi,
1652	. cldh, cldl, cldm, cldt, cldq)
1653	else
1654	CALL nuage (paprs, pplay,
1655	. t_seri, cldliq, cldfra, cldtau, cldemi,
1656	. cldh, cldl, cldm, cldt, cldq)
1657	endif
1658	c
1659	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1660	c
1661	IF (MOD(itaprad,radpas).EQ.0) THEN
1662	DO i = 1, klon
1663	albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
1664	. + falbe(i,is_lic) * pctsrf(i,is_lic)
1665	. + falbe(i,is_ter) * pctsrf(i,is_ter)
1666	. + falbe(i,is_sic) * pctsrf(i,is_sic)
1667	albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
1668	. + falblw(i,is_lic) * pctsrf(i,is_lic)
1669	. + falblw(i,is_ter) * pctsrf(i,is_ter)
1670	. + falblw(i,is_sic) * pctsrf(i,is_sic)
1671	ENDDO
1672	! if (debut) then
1673	! albsol1 = albsol
1674	! albsollw1 = albsollw
1675	! endif
1676	! albsol = albsol1
1677	! albsollw = albsollw1
1678	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
1679	cIM e (dist, rmu0, fract, co2_ppm, solaire,
1680	e (dist, rmu0, fract,
1681	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
1682	e wo,
1683	e cldfra, cldemi, cldtau,
1684	s heat,heat0,cool,cool0,radsol,albpla,
1685	s topsw,toplw,solsw,sollw,
1686	s sollwdown,
1687	cccIMs topsw0,toplw0,solsw0,sollw0)
1688	s topsw0,toplw0,solsw0,sollw0,
1689	s ZFSUP,ZFSDN,ZFSUP0,ZFSDN0)
1690	itaprad = 0
1691	ENDIF
1692	itaprad = itaprad + 1
1693
1694	c
1695	c Ajouter la tendance des rayonnements (tous les pas)
1696	c
1697	DO k = 1, klev
1698	DO i = 1, klon
1699	t_seri(i,k) = t_seri(i,k)
1700	. + (heat(i,k)-cool(i,k)) * dtime/86400.
1701	ENDDO
1702	ENDDO
1703	c
1704	IF (if_ebil.ge.2) THEN
1705	ztit='after rad'
1706	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1707	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1708	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1709	call diagphy(paire,ztit,ip_ebil
1710	e , topsw, toplw, solsw, sollw, zero_v
1711	e , zero_v, zero_v, zero_v, ztsol
1712	e , d_h_vcol, d_qt, d_ec
1713	s , fs_bound, fq_bound )
1714	END IF
1715	c
1716	c
1717	c Calculer l'hydrologie de la surface
1718	c
1719	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
1720	c . agesno, ftsol,fqsol,fsnow, ruis)
1721	c
1722	DO i = 1, klon
1723	zxqsol(i) = 0.0
1724	zxsnow(i) = 0.0
1725	ENDDO
1726	DO nsrf = 1, nbsrf
1727	DO i = 1, klon
1728	zxqsol(i) = zxqsol(i) + fqsol(i,nsrf)*pctsrf(i,nsrf)
1729	zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
1730	ENDDO
1731	ENDDO
1732	c
1733	c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
1734	c
1735	cXXX DO nsrf = 1, nbsrf
1736	cXXX DO i = 1, klon
1737	cXXX IF (pctsrf(i,nsrf).LT.epsfra) THEN
1738	cXXX fqsol(i,nsrf) = zxqsol(i)
1739	cXXX fsnow(i,nsrf) = zxsnow(i)
1740	cXXX ENDIF
1741	cXXX ENDDO
1742	cXXX ENDDO
1743	c
1744	c Calculer le bilan du sol et la derive de temperature (couplage)
1745	c
1746	DO i = 1, klon
1747	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
1748	c a la demande de JLD
1749	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1750	ENDDO
1751	c
1752	cmoddeblott(jan95)
1753	c Appeler le programme de parametrisation de l'orographie
1754	c a l'echelle sous-maille:
1755	c
1756	IF (ok_orodr) THEN
1757	c
1758	c selection des points pour lesquels le shema est actif:
1759	igwd=0
1760	DO i=1,klon
1761	itest(i)=0
1762	c IF ((zstd(i).gt.10.0)) THEN
1763	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1764	itest(i)=1
1765	igwd=igwd+1
1766	idx(igwd)=i
1767	ENDIF
1768	ENDDO
1769	c igwdim=MAX(1,igwd)
1770	c
1771	CALL drag_noro(klon,klev,dtime,paprs,pplay,
1772	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
1773	e igwd,idx,itest,
1774	e t_seri, u_seri, v_seri,
1775	s zulow, zvlow, zustr, zvstr,
1776	s d_t_oro, d_u_oro, d_v_oro)
1777	c
1778	c ajout des tendances
1779	DO k = 1, klev
1780	DO i = 1, klon
1781	t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
1782	u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
1783	v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
1784	ENDDO
1785	ENDDO
1786	c
1787	ENDIF ! fin de test sur ok_orodr
1788	c
1789	IF (ok_orolf) THEN
1790	c
1791	c selection des points pour lesquels le shema est actif:
1792	igwd=0
1793	DO i=1,klon
1794	itest(i)=0
1795	IF ((zpic(i)-zmea(i)).GT.100.) THEN
1796	itest(i)=1
1797	igwd=igwd+1
1798	idx(igwd)=i
1799	ENDIF
1800	ENDDO
1801	c igwdim=MAX(1,igwd)
1802	c
1803	CALL lift_noro(klon,klev,dtime,paprs,pplay,
1804	e rlat,zmea,zstd,zpic,
1805	e itest,
1806	e t_seri, u_seri, v_seri,
1807	s zulow, zvlow, zustr, zvstr,
1808	s d_t_lif, d_u_lif, d_v_lif)
1809	c
1810	c ajout des tendances
1811	DO k = 1, klev
1812	DO i = 1, klon
1813	t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
1814	u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
1815	v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
1816	ENDDO
1817	ENDDO
1818	c
1819	ENDIF ! fin de test sur ok_orolf
1820	c
1821	IF (if_ebil.ge.2) THEN
1822	ztit='after orography'
1823	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1824	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1825	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1826	END IF
1827	c
1828	c
1829	cAA
1830	cAA Installation de l'interface online-offline pour traceurs
1831	cAA
1832	c====================================================================
1833	c Calcul des tendances traceurs
1834	c====================================================================
1835	C Pascale : il faut quand meme apeller phytrac car il gere les sorties
1836	cKE43 des traceurs => il faut donc mettre des flags a .false.
1837	IF (iflag_con.GE.3) THEN
1838	c on ajoute les tendances calculees par KE43
1839	cXXX OM on onhibe la convection sur les traceurs
1840	DO iq=1, nqmax-2 ! Sandrine a -3 ???
1841	cXXX OM on inhibe la convection sur les traceur
1842	cXXX DO k = 1, nlev
1843	cXXX DO i = 1, klon
1844	cXXX tr_seri(i,k,iq) = tr_seri(i,k,iq) + d_tr(i,k,iq)
1845	cXXX ENDDO
1846	cXXX ENDDO
1847	WRITE(iqn,'(i2.2)') iq
1848	CALL minmaxqfi(tr_seri(1,1,iq),0.,1.e33,'couche lim iq='//iqn)
1849	ENDDO
1850	CMAF modif pour garder info du nombre de traceurs auxquels
1851	C la physique s'applique
1852	ELSE
1853	CMAF modif pour garder info du nombre de traceurs auxquels
1854	C la physique s'applique
1855	C
1856	call phytrac (rnpb,
1857	I debut,lafin,
1858	I nqmax-2,
1859	I nlon,nlev,dtime,
1860	I t,paprs,pplay,
1861	I pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1862	I ycoefh,yu1,yv1,ftsol,pctsrf,rlat,
1863	I frac_impa, frac_nucl,
1864	I rlon,presnivs,paire,pphis,
1865	O tr_seri)
1866	ENDIF
1867
1868	IF (offline) THEN
1869
1870	call phystokenc (
1871	I nlon,nlev,pdtphys,rlon,rlat,
1872	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1873	I ycoefh,yu1,yv1,ftsol,pctsrf,
1874	I frac_impa, frac_nucl,
1875	I pphis,paire,dtime,itap)
1876
1877
1878	ENDIF
1879
1880	c
1881	c Calculer le transport de l'eau et de l'energie (diagnostique)
1882	c
1883	CALL transp (paprs,zxtsol,
1884	e t_seri, q_seri, u_seri, v_seri, zphi,
1885	s ve, vq, ue, uq)
1886	c
1887	c Accumuler les variables a stocker dans les fichiers histoire:
1888	c
1889	c
1890	c
1891	c+jld ec_conser
1892	DO k = 1, klev
1893	DO i = 1, klon
1894	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
1895	d_t_ec(i,k)=0.5/ZRCPD
1896	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
1897	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1898	d_t_ec(i,k) = d_t_ec(i,k)/dtime
1899	END DO
1900	END DO
1901	c-jld ec_conser
1902	IF (if_ebil.ge.1) THEN
1903	ztit='after physic'
1904	CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
1905	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1906	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1907	C Comme les tendances de la physique sont ajoute dans la dynamique,
1908	C on devrait avoir que la variation d'entalpie par la dynamique
1909	C est egale a la variation de la physique au pas de temps precedent.
1910	C Donc la somme de ces 2 variations devrait etre nulle.
1911	call diagphy(paire,ztit,ip_ebil
1912	e , topsw, toplw, solsw, sollw, sens
1913	e , evap, rain_fall, snow_fall, ztsol
1914	e , d_h_vcol, d_qt, d_ec
1915	s , fs_bound, fq_bound )
1916	C
1917	d_h_vcol_phy=d_h_vcol
1918	C
1919	END IF
1920	C
1921	cccIM cf. FH
1922	c=======================================================================
1923	c SORTIES
1924	c=======================================================================
1925
1926	c Interpollation sur quelques niveaux de pression
1927	c -----------------------------------------------
1928
1929	call plevel(klon,klev,.true. ,pplay,85000.,u_seri,u850)
1930	call plevel(klon,klev,.false.,pplay,85000.,v_seri,v850)
1931	call plevel(klon,klev,.true. ,pplay,50000.,u_seri,u500)
1932	call plevel(klon,klev,.false.,pplay,50000.,v_seri,v500)
1933	call plevel(klon,klev,.true. ,pplay,20000.,u_seri,u200)
1934	call plevel(klon,klev,.false.,pplay,20000.,v_seri,v200)
1935	call plevel(klon,klev,.true. ,pplay,50000.,zphi,phi500)
1936	call plevel(klon,klev,.true. ,paprs,50000.,omega,w500)
1937
1938	slp(:) = paprs(:,1)exp(pphis(:)/(289.t_seri(:,1)))
1939	c
1940
1941
1942	c=============================================================
1943	c Ecriture des sorties
1944	c=============================================================
1945
1946	#ifdef histhf
1947	#include "write_histhf.h"
1948	#endif
1949
1950	#include "write_histday.h"
1951	#include "write_histmth.h"
1952	#include "write_histins.h"
1953
1954	c=============================================================
1955	c
1956	c Convertir les incrementations en tendances
1957	c
1958	DO k = 1, klev
1959	DO i = 1, klon
1960	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
1961	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
1962	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
1963	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
1964	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
1965	ENDDO
1966	ENDDO
1967	c
1968	IF (nqmax.GE.3) THEN
1969	DO iq = 3, nqmax
1970	DO k = 1, klev
1971	DO i = 1, klon
1972	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
1973	ENDDO
1974	ENDDO
1975	ENDDO
1976	ENDIF
1977	c
1978	c Sauvegarder les valeurs de t et q a la fin de la physique:
1979	c
1980	DO k = 1, klev
1981	DO i = 1, klon
1982	t_ancien(i,k) = t_seri(i,k)
1983	q_ancien(i,k) = q_seri(i,k)
1984	ENDDO
1985	ENDDO
1986	c
1987	c====================================================================
1988	c Si c'est la fin, il faut conserver l'etat de redemarrage
1989	c====================================================================
1990	c
1991	IF (lafin) THEN
1992	itau_phy = itau_phy + itap
1993	ccc IF (ok_oasis) CALL quitcpl
1994	CALL phyredem ("restartphy.nc",dtime,radpas,
1995	. rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsol, fsnow,
1996	. falbe, fevap, rain_fall, snow_fall,
1997	. solsw, sollwdown,dlw,
1998	. radsol,frugs,agesno,
1999	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
2000	. t_ancien, q_ancien, rnebcon, ratqs, clwcon)
2001	ENDIF
2002
2003	RETURN
2004	END
2005	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
2006	IMPLICIT none
2007	c
2008	c Calculer et imprimer l'eau totale. A utiliser pour verifier
2009	c la conservation de l'eau
2010	c
2011	#include "YOMCST.h"
2012	INTEGER klon,klev
2013	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
2014	REAL aire(klon)
2015	REAL qtotal, zx, qcheck
2016	INTEGER i, k
2017	c
2018	zx = 0.0
2019	DO i = 1, klon
2020	zx = zx + aire(i)
2021	ENDDO
2022	qtotal = 0.0
2023	DO k = 1, klev
2024	DO i = 1, klon
2025	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
2026	. *(paprs(i,k)-paprs(i,k+1))/RG
2027	ENDDO
2028	ENDDO
2029	c
2030	qcheck = qtotal/zx
2031	c
2032	RETURN
2033	END
2034	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
2035	IMPLICIT none
2036	c
2037	c Tranformer une variable de la grille physique a
2038	c la grille d'ecriture
2039	c
2040	INTEGER nfield,nlon,iim,jjmp1, jjm
2041	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
2042	c
2043	INTEGER i, n, ig
2044	c
2045	jjm = jjmp1 - 1
2046	DO n = 1, nfield
2047	DO i=1,iim
2048	ecrit(i,n) = fi(1,n)
2049	ecrit(i+jjm*iim,n) = fi(nlon,n)
2050	ENDDO
2051	DO ig = 1, nlon - 2
2052	ecrit(iim+ig,n) = fi(1+ig,n)
2053	ENDDO
2054	ENDDO
2055	RETURN
2056	END
2057

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

Download in other formats:

Original Format