Context Navigation

physiq.F @ 456

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

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

Download in other formats:

Original Format