Context Navigation

physiq.F @ 411

Last change on this file since 411 was 411, checked in by lmdzadmin, 22 years ago

KE vectorise
nouvelles sorties: haute frequence,

rajout diagnostiques pres du sol,
fichiers include pour les sorties,
sorties des flux shortwave au sommet et a la surface

pour l'albedo,

sorties sur niveaux de pression

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

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

Download in other formats:

Original Format