Context Navigation

physiq.F @ 98

Last change on this file since 98 was 97, checked in by slebonnois, 14 years ago
Serie de modifs SL pour homogeneisation des phytitan et phyvenus Ca touche aussi aux liens phy/dyn (surtout a propos de clesphy0), a verifier avec les autres, donc...
File size: 43.3 KB

Line
1	!
2	! $Header: /home/cvsroot/LMDZ4/libf/phylmd/physiq.F,v 1.8 2005/02/24 09:58:18 fairhead Exp $
3	!
4	c
5	SUBROUTINE physiq (nlon,nlev,nqmax,
6	. debut,lafin,rjourvrai,gmtime,pdtphys,
7	. paprs,pplay,ppk,pphi,pphis,presnivs,
8	. u,v,t,qx,
9	. omega,
10	. d_u, d_v, d_t, d_qx, d_ps)
11
12	USE ioipsl
13	USE histcom
14	USE infotrac
15	USE control_mod
16	IMPLICIT none
17	c======================================================================
18	c
19	c Modifications pour la physique de Venus
20	c S. Lebonnois (LMD/CNRS) Septembre 2005
21	c
22	c ---------------------------------------------------------------------
23	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
24	c
25	c Objet: Moniteur general de la physique du modele
26	cAA Modifications quant aux traceurs :
27	cAA - uniformisation des parametrisations ds phytrac
28	cAA - stockage des moyennes des champs necessaires
29	cAA en mode traceur off-line
30	c======================================================================
31	c CLEFS CPP POUR LES IO
32	c =====================
33	c#define histhf
34	#define histday
35	#define histmth
36	#define histins
37	c======================================================================
38	c modif ( P. Le Van , 12/10/98 )
39	c
40	c Arguments:
41	c
42	c nlon----input-I-nombre de points horizontaux
43	c nlev----input-I-nombre de couches verticales
44	c nqmax---input-I-nombre de traceurs
45	c debut---input-L-variable logique indiquant le premier passage
46	c lafin---input-L-variable logique indiquant le dernier passage
47	c rjour---input-R-numero du jour de l'experience
48	c gmtime--input-R-temps universel dans la journee (0 a RDAY s)
49	c pdtphys-input-R-pas d'integration pour la physique (seconde)
50	c paprs---input-R-pression pour chaque inter-couche (en Pa)
51	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
52	c ppk ---input-R-fonction d'Exner au milieu de couche
53	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
54	c pphis---input-R-geopotentiel du sol
55	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
56	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
57	c v-------input-R-vitesse Y (de S a N) en m/s
58	c t-------input-R-temperature (K)
59	c qx------input-R-mass mixing ratio traceurs (kg/kg)
60	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
61	c omega---input-R-vitesse verticale en Pa/s
62	c
63	c d_u-----output-R-tendance physique de "u" (m/s/s)
64	c d_v-----output-R-tendance physique de "v" (m/s/s)
65	c d_t-----output-R-tendance physique de "t" (K/s)
66	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
67	c d_ps----output-R-tendance physique de la pression au sol
68	c======================================================================
69	#include "dimensions.h"
70	integer jjmp1
71	parameter (jjmp1=jjm+1-1/jjm)
72	#include "dimphy.h"
73	#include "dimsoil.h"
74	#include "clesphys.h"
75	#include "temps.h"
76	#include "comgeomphy.h"
77	#include "iniprint.h"
78	#include "raddim.h"
79	#include "timerad.h"
80	#include "logic.h"
81	c======================================================================
82	LOGICAL ok_journe ! sortir le fichier journalier
83	save ok_journe
84	c PARAMETER (ok_journe=.true.)
85	c
86	LOGICAL ok_mensuel ! sortir le fichier mensuel
87	save ok_mensuel
88	c PARAMETER (ok_mensuel=.true.)
89	c
90	LOGICAL ok_instan ! sortir le fichier instantane
91	save ok_instan
92	c PARAMETER (ok_instan=.true.)
93	c
94	c======================================================================
95	c
96	c Variables argument:
97	c
98	INTEGER nlon
99	INTEGER nlev
100	INTEGER nqmax
101	REAL rjourvrai
102	REAL gmtime
103	REAL pdtphys
104	LOGICAL debut, lafin
105	REAL paprs(klon,klev+1)
106	REAL pplay(klon,klev)
107	REAL pphi(klon,klev)
108	REAL pphis(klon)
109	REAL presnivs(klev)
110
111	! ADAPTATION GCM POUR CP(T)
112	REAL ppk(klon,klev)
113
114	REAL u(klon,klev)
115	REAL v(klon,klev)
116	REAL t(klon,klev)
117	REAL qx(klon,klev,nqmax)
118
119	REAL t_ancien(klon,klev)
120	REAL u_ancien(klon,klev)
121	SAVE t_ancien, u_ancien
122	LOGICAL ancien_ok
123	SAVE ancien_ok
124
125	REAL d_u_dyn(klon,klev)
126	REAL d_t_dyn(klon,klev)
127
128	REAL omega(klon,klev)
129
130	REAL d_u(klon,klev)
131	REAL d_v(klon,klev)
132	REAL d_t(klon,klev)
133	REAL d_qx(klon,klev,nqmax)
134	REAL d_ps(klon)
135
136	INTEGER klevp1, klevm1
137	PARAMETER(klevp1=klev+1,klevm1=klev-1)
138	c
139	REAL swnet(klon,klevp1)
140	SAVE swnet
141	c
142	REAL lwnet(klon,klevp1)
143	SAVE lwnet
144	c
145	REAL LWdnTOA(klon)
146	SAVE LWdnTOA
147	c
148	logical ok_hf
149	real ecrit_hf
150	integer nid_hf
151	save ok_hf, ecrit_hf, nid_hf
152
153	#ifdef histhf
154	data ok_hf,ecrit_hf/.true.,0.25/
155	#else
156	data ok_hf/.false./
157	#endif
158
159	c Variables propres a la physique
160	c
161	REAL dtime
162	SAVE dtime ! pas temporel de la physique
163	c
164	INTEGER radpas
165	SAVE radpas ! frequence d'appel rayonnement
166	c
167	REAL radsol(klon)
168	SAVE radsol ! bilan radiatif au sol calcule par code radiatif
169	c
170	REAL rlev(klon,klev+1) ! altitude a chaque niveau (interface inferieure de la couche)
171	SAVE rlev
172	ci
173	INTEGER itap
174	SAVE itap ! compteur pour la physique
175	c
176	REAL ftsol(klon)
177	SAVE ftsol ! temperature du sol
178	c
179	REAL ftsoil(klon,nsoilmx)
180	SAVE ftsoil ! temperature dans le sol
181	c
182	REAL falbe(klon)
183	SAVE falbe ! albedo
184
185	REAL delp(klon,klev) ! epaisseur d'une couche
186
187	CMODDEB FLOTT
188	c
189	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
190	c
191	REAL zmea(klon)
192	SAVE zmea ! orographie moyenne
193	c
194	REAL zstd(klon)
195	SAVE zstd ! deviation standard de l'OESM
196	c
197	REAL zsig(klon)
198	SAVE zsig ! pente de l'OESM
199	c
200	REAL zgam(klon)
201	save zgam ! anisotropie de l'OESM
202	c
203	REAL zthe(klon)
204	SAVE zthe ! orientation de l'OESM
205	c
206	REAL zpic(klon)
207	SAVE zpic ! Maximum de l'OESM
208	c
209	REAL zval(klon)
210	SAVE zval ! Minimum de l'OESM
211	c
212	REAL rugoro(klon)
213	SAVE rugoro ! longueur de rugosite de l'OESM
214
215	INTEGER igwd,idx(klon),itest(klon)
216	c
217	c Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro
218
219	REAL zulow(klon),zvlow(klon)
220	REAL zustrdr(klon), zvstrdr(klon)
221	REAL zustrli(klon), zvstrli(klon)
222	REAL zustrhi(klon), zvstrhi(klon)
223
224	c Pour calcul GW drag oro et nonoro: CALCUL de N2:
225	real zdzlev(klon,klev)
226	real ztlev(klon,klev),zpklev(klon,klev)
227	real ztetalay(klon,klev),ztetalev(klon,klev)
228	real zdtetalev(klon,klev)
229	real zn2(klon,klev) ! BV^2 at plev
230
231	c Pour les bilans de moment angulaire,
232	integer bilansmc
233	c Pour le transport de ballons
234	integer ballons
235	c j'ai aussi besoin
236	c du stress de couche limite a la surface:
237
238	REAL zustrcl(klon),zvstrcl(klon)
239
240	c et du stress total c de la physique:
241
242	REAL zustrph(klon),zvstrph(klon)
243	c
244	REAL zuthe(klon),zvthe(klon)
245	SAVE zuthe
246	SAVE zvthe
247
248	c Variables locales:
249	c
250	REAL cdragh(klon) ! drag coefficient pour T and Q
251	REAL cdragm(klon) ! drag coefficient pour vent
252	c
253	cAA Pour TRACEURS
254	cAA
255	REAL source(klon,nqmax)
256	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
257	REAL yu1(klon) ! vents dans la premiere couche U
258	REAL yv1(klon) ! vents dans la premiere couche V
259
260	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
261	REAL dlw(klon) ! derivee infra rouge
262	cym
263	SAVE dlw
264	cym
265	REAL fder(klon) ! Derive de flux (sensible et latente)
266	save fder
267	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
268	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
269	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
270	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
271	c
272
273	c======================================================================
274	c
275	c Declaration des procedures appelees
276	c
277	EXTERNAL ajsec ! ajustement sec
278	EXTERNAL clmain ! couche limite
279	EXTERNAL hgardfou ! verifier les temperatures
280	c EXTERNAL orbite ! calculer l'orbite
281	EXTERNAL phyetat0 ! lire l'etat initial de la physique
282	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
283	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
284	EXTERNAL suphec ! initialiser certaines constantes
285	EXTERNAL transp ! transport total de l'eau et de l'energie
286	INTEGER lnblnk1
287	EXTERNAL lnblnk1 !enleve les blancs a la fin d'une variable de type
288	!caracter
289	EXTERNAL abort_gcm
290	EXTERNAL printflag
291	EXTERNAL zenang
292	EXTERNAL diagetpq
293	EXTERNAL conf_phys
294	EXTERNAL diagphy
295	EXTERNAL mucorr
296	EXTERNAL phytrac
297	c
298	c Variables locales
299	c
300	CXXX PB
301	REAL fluxt(klon,klev) ! flux turbulent de chaleur
302	REAL fluxu(klon,klev) ! flux turbulent de vitesse u
303	REAL fluxv(klon,klev) ! flux turbulent de vitesse v
304	c
305	REAL flux_dyn(klon,klev) ! flux de chaleur produit par la dynamique
306	REAL flux_ajs(klon,klev) ! flux de chaleur ajustement sec
307	REAL flux_ec(klon,klev) ! flux de chaleur Ec
308	c
309	REAL heat(klon,klev) ! chauffage solaire
310	REAL cool(klon,klev) ! refroidissement infrarouge
311	REAL dtrad(klon,klev) ! K s-1
312	REAL tmpout(klon,klev) ! K s-1
313	REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)
314	real sollwdown(klon) ! downward LW flux at surface
315
316	c Le rayonnement n'est pas calcule tous les pas, il faut donc
317	c sauvegarder les sorties du rayonnement
318	SAVE heat,cool,dtrad,topsw,toplw,solsw,sollw,sollwdown
319	c
320	INTEGER itaprad
321	SAVE itaprad
322	c
323	REAL dist, rmu0(klon), fract(klon)
324	REAL zdtime, zlongi
325	c
326	INTEGER i, k, iq, ig, j, ll
327	c
328	REAL zphi(klon,klev)
329	c
330	c Variables du changement
331	c
332	c ajs: ajustement sec
333	c vdf: couche limite (Vertical DiFfusion)
334	REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
335	REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
336	c
337	REAL d_ts(klon)
338	c
339	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
340	REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
341	c
342	CMOD LOTT: Tendances Orography Sous-maille
343	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
344	REAL d_t_oro(klon,klev)
345	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
346	REAL d_t_lif(klon,klev)
347	C Tendances Ondes de G non oro (runs strato).
348	REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
349	REAL d_t_hin(klon,klev)
350
351	c
352	c Variables liees a l'ecriture de la bande histoire physique
353	c
354	INTEGER ecrit_mth
355	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
356	c
357	INTEGER ecrit_day
358	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
359	c
360	INTEGER ecrit_ins
361	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
362	c
363	integer itau_w ! pas de temps ecriture = itap + itau_phy
364
365	c Variables locales pour effectuer les appels en serie
366	c
367	REAL t_seri(klon,klev)
368	REAL u_seri(klon,klev), v_seri(klon,klev)
369	c
370	REAL tr_seri(klon,klev,nqmax)
371	REAL d_tr(klon,klev,nqmax)
372
373	INTEGER length
374	PARAMETER ( length = 100 )
375	REAL tabcntr0( length )
376	c
377	c pour ioipsl
378	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
379	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
380	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
381	REAL zx_tmp_2d(iim,jjmp1),zx_tmp_3d(iim,jjmp1,klev)
382	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
383
384	INTEGER nid_day, nid_mth, nid_ins
385	SAVE nid_day, nid_mth, nid_ins
386	c
387	INTEGER nhori, nvert, idayref
388	REAL zsto, zout, zsto1, zsto2, zero
389	parameter (zero=0.0e0)
390	real zjulian
391	save zjulian
392
393	CHARACTER*2 str2
394	character*20 modname
395	character*80 abort_message
396	logical ok_sync
397
398	character*30 nom_fichier
399	character*10 varname
400	character*40 vartitle
401	character*20 varunits
402	C Variables liees au bilan d'energie et d'enthalpi
403	REAL ztsol(klon)
404	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
405	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
406	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
407	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
408	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
409	REAL d_h_vcol_phy
410	REAL fs_bound, fq_bound
411	SAVE d_h_vcol_phy
412	REAL zero_v(klon),zero_v2(klon,klev)
413	CHARACTER*15 ztit
414	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
415	SAVE ip_ebil
416	DATA ip_ebil/2/
417	INTEGER if_ebil ! level for energy conserv. dignostics
418	SAVE if_ebil
419	c+jld ec_conser
420	REAL d_t_ec(klon,klev) ! tendance du a la conversion Ec -> E thermique
421	c-jld ec_conser
422
423	c TEST VENUS...
424	REAL mang(klon,klev) ! moment cinetique
425	REAL mangtot ! moment cinetique total
426
427	c Declaration des constantes et des fonctions thermodynamiques
428	c
429	#include "YOMCST.h"
430
431	c======================================================================
432	c INITIALISATIONS
433	c================
434
435	modname = 'physiq'
436	ok_sync=.TRUE.
437
438	bilansmc = 1
439	ballons = 0
440
441	IF (if_ebil.ge.1) THEN
442	DO i=1,klon
443	zero_v(i)=0.
444	END DO
445	DO i=1,klon
446	DO j=1,klev
447	zero_v2(i,j)=0.
448	END DO
449	END DO
450	END IF
451
452	c PREMIER APPEL SEULEMENT
453	c========================
454	IF (debut) THEN
455
456	CALL suphec ! initialiser constantes et parametres phys.
457
458	IF (if_ebil.ge.1) d_h_vcol_phy=0.
459	c
460	c appel a la lecture du physiq.def
461	c
462	call conf_phys(ok_journe, ok_mensuel,
463	. ok_instan,
464	. if_ebil)
465
466	c
467	c Initialiser les compteurs:
468	c
469	itap = 0
470	itaprad = 0
471	c
472	c Lecture startphy.nc :
473	c
474	c REMETTRE TOUS LES PARAMETRES POUR OROGW...
475	CALL phyetat0 ("startphy.nc",dtime,
476	. rlatd,rlond,ftsol,ftsoil,
477	. falbe, solsw, sollw,
478	. dlw,radsol,
479	. zmea,zstd,zsig,zgam,zthe,zpic,zval,
480	. tabcntr0,
481	. t_ancien, ancien_ok)
482
483	radpas = NINT( RDAY/pdtphys/nbapp_rad)
484
485	CALL printflag( tabcntr0,radpas,ok_journe,ok_instan )
486	c
487	c---------
488	c FLOTT
489	IF (ok_orodr) THEN
490	DO i=1,klon
491	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
492	ENDDO
493	CALL SUGWD(klon,klev,paprs,pplay)
494	DO i=1,klon
495	zuthe(i)=0.
496	zvthe(i)=0.
497	if(zstd(i).gt.10.)then
498	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
499	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
500	endif
501	ENDDO
502	ENDIF
503
504	if (bilansmc.eq.1) then
505	C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
506	C DU BILAN DE MOMENT ANGULAIRE.
507	open(27,file='aaam_bud.out',form='formatted')
508	open(28,file='fields_2d.out',form='formatted')
509	write(,)'Ouverture de aaam_bud.out (FL Vous parle)'
510	write(,)'Ouverture de fields_2d.out (FL Vous parle)'
511	endif !bilansmc
512
513	c--------------SLEBONNOIS
514	C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
515	C DES BALLONS
516	if (ballons.eq.1) then
517	open(30,file='ballons-lat.out',form='formatted')
518	open(31,file='ballons-lon.out',form='formatted')
519	open(32,file='ballons-u.out',form='formatted')
520	open(33,file='ballons-v.out',form='formatted')
521	open(34,file='ballons-alt.out',form='formatted')
522	write(,)'Ouverture des ballons*.out'
523	endif !ballons
524	c-------------
525
526	c---------
527	C TRACEURS
528	C source dans couche limite
529	source = 0.0 ! pas de source, pour l'instant
530	C
531	c---------
532	c
533	c Verifications:
534	c
535	IF (ABS(dtime-pdtphys).GT.0.001) THEN
536	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
537	. pdtphys
538	c abort_message='Pas physique n est pas correct '
539	c call abort_gcm(modname,abort_message,1)
540	dtime=pdtphys
541	ENDIF
542	IF (nlon .NE. klon) THEN
543	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
544	. klon
545	abort_message='nlon et klon ne sont pas coherents'
546	call abort_gcm(modname,abort_message,1)
547	ENDIF
548	IF (nlev .NE. klev) THEN
549	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
550	. klev
551	abort_message='nlev et klev ne sont pas coherents'
552	call abort_gcm(modname,abort_message,1)
553	ENDIF
554	c
555	IF (dtimeFLOAT(radpas).GT.(RDAY0.25).AND.cycle_diurne)
556	$ THEN
557	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
558	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
559	abort_message='Nbre d appels au rayonnement insuffisant'
560	call abort_gcm(modname,abort_message,1)
561	ENDIF
562	c
563	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
564	. iflag_ajs
565	c
566	ecrit_mth = NINT(RDAY/dtime*ecritphy) ! tous les ecritphy jours
567	IF (ok_mensuel) THEN
568	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
569	. ecrit_mth
570	ENDIF
571
572	ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours
573	IF (ok_journe) THEN
574	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
575	. ecrit_day
576	ENDIF
577
578	ecrit_ins = NINT(RDAY/dtime/*ecritphy) ! Fraction de jour reglable
579	IF (ok_instan) THEN
580	WRITE(lunout,*)'La frequence de sortie instant. est de ',
581	. ecrit_ins
582	ENDIF
583
584	c Initialisation des sorties
585	c===========================
586
587	#ifdef CPP_IOIPSL
588
589	#ifdef histhf
590	#include "ini_histhf.h"
591	#endif
592
593	#ifdef histday
594	#include "ini_histday.h"
595	#endif
596
597	#ifdef histmth
598	#include "ini_histmth.h"
599	#endif
600
601	#ifdef histins
602	#include "ini_histins.h"
603	#endif
604
605	#endif
606
607	c
608	c Initialiser les valeurs de u pour calculs tendances
609	c (pour T, c'est fait dans phyetat0)
610	c
611	DO k = 1, klev
612	DO i = 1, klon
613	u_ancien(i,k) = u(i,k)
614	ENDDO
615	ENDDO
616
617
618	ENDIF ! debut
619	c====================================================================
620	c======================================================================
621
622	c Mettre a zero des variables de sortie (pour securite)
623	c
624	DO i = 1, klon
625	d_ps(i) = 0.0
626	ENDDO
627	DO k = 1, klev
628	DO i = 1, klon
629	d_t(i,k) = 0.0
630	d_u(i,k) = 0.0
631	d_v(i,k) = 0.0
632	ENDDO
633	ENDDO
634	DO iq = 1, nqmax
635	DO k = 1, klev
636	DO i = 1, klon
637	d_qx(i,k,iq) = 0.0
638	ENDDO
639	ENDDO
640	ENDDO
641	c
642	c Ne pas affecter les valeurs entrees de u, v, h, et q
643	c
644	DO k = 1, klev
645	DO i = 1, klon
646	t_seri(i,k) = t(i,k)
647	u_seri(i,k) = u(i,k)
648	v_seri(i,k) = v(i,k)
649	ENDDO
650	ENDDO
651	DO iq = 1, nqmax
652	DO k = 1, klev
653	DO i = 1, klon
654	tr_seri(i,k,iq) = qx(i,k,iq)
655	ENDDO
656	ENDDO
657	ENDDO
658	C
659	DO i = 1, klon
660	ztsol(i) = ftsol(i)
661	ENDDO
662	C
663	IF (if_ebil.ge.1) THEN
664	ztit='after dynamic'
665	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
666	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
667	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
668	C Comme les tendances de la physique sont ajoute dans la dynamique,
669	C on devrait avoir que la variation d'entalpie par la dynamique
670	C est egale a la variation de la physique au pas de temps precedent.
671	C Donc la somme de ces 2 variations devrait etre nulle.
672	call diagphy(airephy,ztit,ip_ebil
673	e , zero_v, zero_v, zero_v, zero_v, zero_v
674	e , zero_v, zero_v, zero_v, ztsol
675	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
676	s , fs_bound, fq_bound )
677	END IF
678
679	c====================================================================
680	c Diagnostiquer la tendance dynamique
681	c
682	IF (ancien_ok) THEN
683	DO k = 1, klev
684	DO i = 1, klon
685	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
686	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
687	ENDDO
688	ENDDO
689
690	! ADAPTATION GCM POUR CP(T)
691	do i=1,klon
692	flux_dyn(i,1) = 0.0
693	do j=2,klev
694	flux_dyn(i,j) = flux_dyn(i,j-1)
695	. +cpdet(t_seri(i,j-1))/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
696	enddo
697	enddo
698
699	ELSE
700	DO k = 1, klev
701	DO i = 1, klon
702	d_u_dyn(i,k) = 0.0
703	d_t_dyn(i,k) = 0.0
704	ENDDO
705	ENDDO
706	ancien_ok = .TRUE.
707	ENDIF
708	c====================================================================
709	c
710	c Ajouter le geopotentiel du sol:
711	c
712	DO k = 1, klev
713	DO i = 1, klon
714	zphi(i,k) = pphi(i,k) + pphis(i)
715	ENDDO
716	ENDDO
717	c====================================================================
718	c
719	c Verifier les temperatures
720	c
721	CALL hgardfou(t_seri,ftsol,'debutphy')
722	c====================================================================
723	c
724	c Incrementer le compteur de la physique
725	c
726	itap = itap + 1
727
728	c====================================================================
729	c Orbite et eclairement
730	c====================================================================
731
732	c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
733	c donc pas de variations de Ls, ni de dist.
734	c La seule chose qui compte, c'est la rotation de la planete devant
735	c le Soleil...
736
737	zlongi = 0.0
738	dist = 0.72 ! en UA
739
740	c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
741	c a sa valeur, et prendre en compte leur evolution,
742	c il faudra refaire un orbite.F...
743	c CALL orbite(zlongi,dist)
744
745	IF (cycle_diurne) THEN
746	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
747	CALL zenang(zlongi,gmtime,zdtime,rlatd,rlond,rmu0,fract)
748	ELSE
749	call mucorr(klon,zlongi,rlatd,rmu0,fract)
750	ENDIF
751
752	c====================================================================
753	c Calcul des tendances traceurs
754	c====================================================================
755
756	if (iflag_trac.eq.1) then
757	call phytrac (
758	I itap, gmtime,
759	I debut,lafin,
760	I nqmax,
761	I nlon,nlev,dtime,
762	I u,v,t,paprs,pplay,
763	I rlatd,
764	I rlond,presnivs,pphis,pphi,
765	I falbe,
766	O tr_seri)
767	endif
768
769	c
770	c====================================================================
771	c Appeler la diffusion verticale (programme de couche limite)
772	c====================================================================
773
774	c-------------------------------
775	c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force
776	c l'equilibre radiatif du sol
777	if (1.eq.0) then
778	if (debut) then
779	print*,"ATTENTION, CLMAIN SHUNTEE..."
780	endif
781
782	DO i = 1, klon
783	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
784	fder(i) = 0.0e0
785	dlw(i) = 0.0e0
786	ENDDO
787
788	c Incrementer la temperature du sol
789	c
790	DO i = 1, klon
791	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
792	ftsol(i) = ftsol(i) + d_ts(i)
793	do j=1,nsoilmx
794	ftsoil(i,j)=ftsol(i)
795	enddo
796	ENDDO
797
798	c-------------------------------
799	else
800	c-------------------------------
801
802	fder = dlw
803
804	c print*,"radsol avant clmain=",radsol(klon/2)
805	c print*,"solsw avant clmain=",solsw(klon/2)
806	c print*,"sollw avant clmain=",sollw(klon/2)
807
808	! ADAPTATION GCM POUR CP(T)
809	CALL clmain(dtime,itap,
810	e t_seri,u_seri,v_seri,
811	e rmu0,
812	e ftsol,
813	$ ftsoil,
814	$ paprs,pplay,ppk,radsol,falbe,
815	e solsw, sollw, sollwdown, fder,
816	e rlond, rlatd, cuphy, cvphy,
817	e debut, lafin,
818	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
819	s fluxt,fluxu,fluxv,cdragh,cdragm,
820	s dsens,
821	s ycoefh,yu1,yv1)
822
823	c print*,"radsol apres clmain=",radsol(klon/2)
824	c print*,"solsw apres clmain=",solsw(klon/2)
825	c print*,"sollw apres clmain=",sollw(klon/2)
826
827	CXXX Incrementation des flux
828	DO i = 1, klon
829	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
830	fder(i) = dlw(i) + dsens(i)
831	ENDDO
832	CXXX
833
834	DO k = 1, klev
835	DO i = 1, klon
836	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
837	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
838	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
839	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
840	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
841	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
842	ENDDO
843	ENDDO
844	c
845	c print,"d_t_vdf1=",d_t_vdf(1,:)dtime
846	c print,"d_t_vdf2=",d_t_vdf(klon/2,:)dtime
847	c print,"d_t_vdf3=",d_t_vdf(klon,:)dtime
848	c print,"d_u_vdf=",d_u_vdf(klon/2,:)dtime
849	c print,"d_v_vdf=",d_v_vdf(klon/2,:)dtime
850
851	C TRACEURS
852
853	if (iflag_trac.eq.1) then
854	DO iq=1, nqmax
855	CALL cltrac(dtime,ycoefh,t_seri,
856	s tr_seri(1,1,iq),source,
857	e paprs, pplay,delp,
858	s d_tr_vdf(1,1,iq))
859	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
860	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
861	ENDDO
862	endif
863
864	IF (if_ebil.ge.2) THEN
865	ztit='after clmain'
866	CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
867	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
868	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
869	call diagphy(airephy,ztit,ip_ebil
870	e , zero_v, zero_v, zero_v, zero_v, sens
871	e , zero_v, zero_v, zero_v, ztsol
872	e , d_h_vcol, d_qt, d_ec
873	s , fs_bound, fq_bound )
874	END IF
875	C
876	c
877	c Incrementer la temperature du sol
878	c
879	c print*,'Tsol avant clmain:',ftsol(1)
880	DO i = 1, klon
881	ftsol(i) = ftsol(i) + d_ts(i)
882	ENDDO
883	c print*,'Tsol apres clmain:',ftsol(1)
884
885	c Calculer la derive du flux infrarouge
886	c
887	DO i = 1, klon
888	dlw(i) = - 4.0RSIGMAftsol(i)**3
889	ENDDO
890
891	c-------------------------------
892	endif ! fin du VENUS TEST
893
894	c
895	c Appeler l'ajustement sec
896	c
897	c===================================================================
898	c Convection seche
899	c===================================================================
900	c
901	d_t_ajs(:,:)=0.
902	d_u_ajs(:,:)=0.
903	d_v_ajs(:,:)=0.
904	d_tr_ajs(:,:,:)=0.
905	c
906	IF(prt_level>9)WRITE(lunout,*)
907	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
908	s ,iflag_ajs
909
910	if(iflag_ajs.eq.0) then
911	c Rien
912	c ====
913	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
914
915	else if(iflag_ajs.eq.1) then
916
917	c Ajustement sec
918	c ==============
919	IF(prt_level>9)WRITE(lunout,*)'ajsec'
920
921	! ADAPTATION GCM POUR CP(T)
922	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
923	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
924
925	! ADAPTATION GCM POUR CP(T)
926	do i=1,klon
927	flux_ajs(i,1) = 0.0
928	do j=2,klev
929	flux_ajs(i,j) = flux_ajs(i,j-1)
930	. + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
931	. *(paprs(i,j-1)-paprs(i,j))
932	enddo
933	enddo
934
935	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
936	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
937	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
938	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
939	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
940	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
941	if (iflag_trac.eq.1) then
942	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
943	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
944	endif
945
946	c print,"d_t_ajs1=",d_t_ajs(1,:)dtime
947	c print,"d_t_ajs2=",d_t_ajs(klon/2,:)dtime
948	c print,"d_t_ajs3=",d_t_ajs(klon,:)dtime
949	c print,"d_u_ajs=",d_u_ajs(klon/2,:)dtime
950	c print,"d_v_ajs=",d_v_ajs(klon/2,:)dtime
951
952	endif
953	c
954	IF (if_ebil.ge.2) THEN
955	ztit='after dry_adjust'
956	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
957	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
958	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
959	call diagphy(airephy,ztit,ip_ebil
960	e , zero_v, zero_v, zero_v, zero_v, sens
961	e , zero_v, zero_v, zero_v, ztsol
962	e , d_h_vcol, d_qt, d_ec
963	s , fs_bound, fq_bound )
964	END IF
965
966	c====================================================================
967	c RAYONNEMENT
968	c====================================================================
969
970	IF (MOD(itaprad,radpas).EQ.0) THEN
971	c print*,'RAYONNEMENT ',
972	c $ ' (itaprad=',itaprad,'/radpas=',radpas,')'
973
974	dtimerad = dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
975	c PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
976
977	c print*,"radsol avant radlwsw=",radsol(klon/2)
978	c print*,"solsw avant radlwsw=",solsw(klon/2)
979	c print*,"sollw avant radlwsw=",sollw(klon/2)
980	c print*,"avant radlwsw"
981	CALL radlwsw
982	e (dist, rmu0, fract,
983	e paprs, pplay,ftsol, t_seri,
984	s heat,cool,radsol,
985	s topsw,toplw,solsw,sollw,
986	s sollwdown,
987	s lwnet, swnet)
988	c print*,"apres radlwsw"
989
990	c print*,"radsol apres radlwsw=",radsol(klon/2)
991	c print*,"solsw apres radlwsw=",solsw(klon/2)
992	c print*,"sollw apres radlwsw=",sollw(klon/2)
993	itaprad = 0
994	DO k = 1, klev
995	DO i = 1, klon
996	dtrad(i,k) = heat(i,k)-cool(i,k)
997	dtrad(i,k) = dtrad(i,k)/RDAY !K/s
998	ENDDO
999	ENDDO
1000	c print,"dtrad1=",dtrad(1,:)dtime
1001	c print,"dtrad2=",dtrad(klon/2,:)dtime
1002	c print,"dtrad3=",dtrad(klon,:)dtime
1003
1004	ENDIF
1005	itaprad = itaprad + 1
1006	c====================================================================
1007	c
1008	c Ajouter la tendance des rayonnements (tous les pas)
1009	c
1010	DO k = 1, klev
1011	DO i = 1, klon
1012	t_seri(i,k) = t_seri(i,k) + dtrad(i,k) * dtime
1013	ENDDO
1014	ENDDO
1015
1016	IF (if_ebil.ge.2) THEN
1017	ztit='after rad'
1018	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1019	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1020	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1021	call diagphy(airephy,ztit,ip_ebil
1022	e , topsw, toplw, solsw, sollw, zero_v
1023	e , zero_v, zero_v, zero_v, ztsol
1024	e , d_h_vcol, d_qt, d_ec
1025	s , fs_bound, fq_bound )
1026	END IF
1027	c
1028
1029	c====================================================================
1030	c Calcul des gravity waves FLOTT
1031	c====================================================================
1032	c
1033	if (ok_orodr.or.ok_gw_nonoro) then
1034	c CALCUL DE N2
1035	do i=1,klon
1036	do k=2,klev
1037	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1038	zpklev(i,k) = sqrt(ppk(i,k)*ppk(i,k-1))
1039	enddo
1040	enddo
1041	call t2tpot(klon*klev,ztlev, ztetalev,zpklev)
1042	call t2tpot(klon*klev,t_seri,ztetalay,ppk)
1043	do i=1,klon
1044	do k=2,klev
1045	zdtetalev(i,k) = ztetalay(i,k)-ztetalay(i,k-1)
1046	zdzlev(i,k) = (zphi(i,k)-zphi(i,k-1))/RG
1047	zn2(i,k) = RGzdtetalev(i,k)/(ztetalev(i,k)zdzlev(i,k))
1048	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
1049	enddo
1050	enddo
1051
1052	endif
1053
1054	c ----------------------------ORODRAG
1055	IF (ok_orodr) THEN
1056	c
1057	c selection des points pour lesquels le shema est actif:
1058	igwd=0
1059	DO i=1,klon
1060	itest(i)=0
1061	c IF ((zstd(i).gt.10.0)) THEN
1062	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1063	itest(i)=1
1064	igwd=igwd+1
1065	idx(igwd)=i
1066	ENDIF
1067	ENDDO
1068	c igwdim=MAX(1,igwd)
1069	c
1070	c A ADAPTER POUR VENUS!!!
1071	CALL drag_noro(klon,klev,dtime,paprs,pplay,zphi,zn2,
1072	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
1073	e igwd,idx,itest,
1074	e t_seri, u_seri, v_seri,
1075	s zulow, zvlow, zustrdr, zvstrdr,
1076	s d_t_oro, d_u_oro, d_v_oro)
1077
1078	c ajout des tendances
1079	t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
1080	d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
1081	u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
1082	d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
1083	v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
1084	d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
1085	c
1086	ELSE
1087	d_t_oro = 0.
1088	d_u_oro = 0.
1089	d_v_oro = 0.
1090	zustrdr = 0.
1091	zvstrdr = 0.
1092	c
1093	ENDIF ! fin de test sur ok_orodr
1094	c
1095	c ----------------------------OROLIFT
1096	IF (ok_orolf) THEN
1097	c
1098	c selection des points pour lesquels le shema est actif:
1099	igwd=0
1100	DO i=1,klon
1101	itest(i)=0
1102	IF ((zpic(i)-zmea(i)).GT.100.) THEN
1103	itest(i)=1
1104	igwd=igwd+1
1105	idx(igwd)=i
1106	ENDIF
1107	ENDDO
1108	c igwdim=MAX(1,igwd)
1109	c
1110	c A ADAPTER POUR VENUS!!!
1111	c CALL lift_noro(klon,klev,dtime,paprs,pplay,
1112	c e rlatd,zmea,zstd,zpic,zgam,zthe,zpic,zval,
1113	c e igwd,idx,itest,
1114	c e t_seri, u_seri, v_seri,
1115	c s zulow, zvlow, zustrli, zvstrli,
1116	c s d_t_lif, d_u_lif, d_v_lif )
1117
1118	c
1119	c ajout des tendances
1120	t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
1121	d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
1122	u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
1123	d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
1124	v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
1125	d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
1126	c
1127	ELSE
1128	d_t_lif = 0.
1129	d_u_lif = 0.
1130	d_v_lif = 0.
1131	zustrli = 0.
1132	zvstrli = 0.
1133	c
1134	ENDIF ! fin de test sur ok_orolf
1135
1136	c ---------------------------- NON-ORO GRAVITY WAVES
1137	IF(ok_gw_nonoro) then
1138
1139	call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
1140	e t_seri, u_seri, v_seri,
1141	o zustrhi,zvstrhi,
1142	o d_t_hin, d_u_hin, d_v_hin)
1143
1144	c ajout des tendances
1145
1146	t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
1147	d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
1148	u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
1149	d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
1150	v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
1151	d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
1152
1153	ELSE
1154	d_t_hin = 0.
1155	d_u_hin = 0.
1156	d_v_hin = 0.
1157	zustrhi = 0.
1158	zvstrhi = 0.
1159
1160	ENDIF ! fin de test sur ok_gw_nonoro
1161
1162	c====================================================================
1163	c Transport de ballons
1164	c====================================================================
1165	if (ballons.eq.1) then
1166	CALL ballon(30,pdtphys,rjourvrai,gmtime,rlatd,rlond,
1167	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
1168	C t,pplay,u,v,zphi) ! alt above planet average radius
1169	endif !ballons
1170
1171	c====================================================================
1172	c Bilan de mmt angulaire
1173	c====================================================================
1174	if (bilansmc.eq.1) then
1175	CMODDEB FLOTT
1176	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
1177	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
1178
1179	DO i = 1, klon
1180	zustrph(i)=0.
1181	zvstrph(i)=0.
1182	zustrcl(i)=0.
1183	zvstrcl(i)=0.
1184	ENDDO
1185	DO k = 1, klev
1186	DO i = 1, klon
1187	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
1188	c (paprs(i,k)-paprs(i,k+1))/rg
1189	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
1190	c (paprs(i,k)-paprs(i,k+1))/rg
1191	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
1192	c (paprs(i,k)-paprs(i,k+1))/rg
1193	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
1194	c (paprs(i,k)-paprs(i,k+1))/rg
1195	ENDDO
1196	ENDDO
1197
1198	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
1199	C ra,rg,romega,
1200	C rlatd,rlond,pphis,
1201	C zustrdr,zustrli,zustrcl,
1202	C zvstrdr,zvstrli,zvstrcl,
1203	C paprs,u,v)
1204
1205	CCMODFIN FLOTT
1206	endif !bilansmc
1207
1208	c====================================================================
1209	c====================================================================
1210	c Calculer le transport de l'eau et de l'energie (diagnostique)
1211	c
1212	c A REVOIR POUR VENUS...
1213	c
1214	c CALL transp (paprs,ftsol,
1215	c e t_seri, q_seri, u_seri, v_seri, zphi,
1216	c s ve, vq, ue, uq)
1217	c
1218	c====================================================================
1219	c+jld ec_conser
1220	DO k = 1, klev
1221	DO i = 1, klon
1222	d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
1223	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
1224	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1225	d_t_ec(i,k) = d_t_ec(i,k)/dtime
1226	END DO
1227	END DO
1228	do i=1,klon
1229	flux_ec(i,1) = 0.0
1230	do j=2,klev
1231	flux_ec(i,j) = flux_ec(i,j-1)
1232	. +cpdet(t_seri(i,j-1))/RGd_t_ec(i,j-1)(paprs(i,j-1)-paprs(i,j))
1233	enddo
1234	enddo
1235
1236	c-jld ec_conser
1237	c====================================================================
1238	IF (if_ebil.ge.1) THEN
1239	ztit='after physic'
1240	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
1241	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1242	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1243	C Comme les tendances de la physique sont ajoute dans la dynamique,
1244	C on devrait avoir que la variation d'entalpie par la dynamique
1245	C est egale a la variation de la physique au pas de temps precedent.
1246	C Donc la somme de ces 2 variations devrait etre nulle.
1247	call diagphy(airephy,ztit,ip_ebil
1248	e , topsw, toplw, solsw, sollw, sens
1249	e , zero_v, zero_v, zero_v, ztsol
1250	e , d_h_vcol, d_qt, d_ec
1251	s , fs_bound, fq_bound )
1252	C
1253	d_h_vcol_phy=d_h_vcol
1254	C
1255	END IF
1256	C
1257	c=======================================================================
1258	c SORTIES
1259	c=======================================================================
1260
1261	c Convertir les incrementations en tendances
1262	c
1263	DO k = 1, klev
1264	DO i = 1, klon
1265	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
1266	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
1267	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
1268	ENDDO
1269	ENDDO
1270	c
1271	DO iq = 1, nqmax
1272	DO k = 1, klev
1273	DO i = 1, klon
1274	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
1275	ENDDO
1276	ENDDO
1277	ENDDO
1278
1279	c------------------------
1280	c Calcul moment cinetique
1281	c------------------------
1282	c TEST VENUS...
1283	c mangtot = 0.0
1284	c DO k = 1, klev
1285	c DO i = 1, klon
1286	c mang(i,k) = RAcos(rlatd(i)RPI/180.)
1287	c . (u_seri(i,k)+RAcos(rlatd(i)RPI/180.)ROMEGA)
1288	c . airephy(i)(paprs(i,k)-paprs(i,k+1))/RG
1289	c mangtot=mangtot+mang(i,k)
1290	c ENDDO
1291	c ENDDO
1292	c print*,"Moment cinetique total = ",mangtot
1293	c
1294	c------------------------
1295	c
1296	c Sauvegarder les valeurs de t et u a la fin de la physique:
1297	c
1298	DO k = 1, klev
1299	DO i = 1, klon
1300	u_ancien(i,k) = u_seri(i,k)
1301	t_ancien(i,k) = t_seri(i,k)
1302	ENDDO
1303	ENDDO
1304	c
1305	c=============================================================
1306	c Ecriture des sorties
1307	c=============================================================
1308
1309	#ifdef CPP_IOIPSL
1310
1311	#ifdef histhf
1312	#include "write_histhf.h"
1313	#endif
1314
1315	#ifdef histday
1316	#include "write_histday.h"
1317	#endif
1318
1319	#ifdef histmth
1320	#include "write_histmth.h"
1321	#endif
1322
1323	#ifdef histins
1324	#include "write_histins.h"
1325	#endif
1326
1327	#endif
1328
1329	c ---------- Sorties testphys1d -------------
1330
1331	if (klon.eq.1) then
1332	call writeg1d(klon,klev,t_seri,"Temp","Temperature")
1333	call writeg1d(klon,1,ftsol,"tsurf","Surface Temp")
1334	DO k = 1, klev
1335	DO i = 1, klon
1336	c tmpout(i,k) = real(heat(i,k))/RDAY
1337	tmpout(i,k) = heat(i,k)/RDAY
1338	ENDDO
1339	ENDDO
1340	call writeg1d(klon,klev,tmpout,
1341	. "heat","Solar heating")
1342	DO k = 1, klev
1343	DO i = 1, klon
1344	c tmpout(i,k) = real(dtrad(i,k))
1345	tmpout(i,k) = dtrad(i,k)
1346	ENDDO
1347	ENDDO
1348	call writeg1d(klon,klev,tmpout,
1349	. "dtrad","IR cooling")
1350	call writeg1d(klon,klev,lwnet,"lwnet","Net LW flux")
1351	call writeg1d(klon,klev,swnet,"swnet","Net SW flux")
1352	call writeg1d(klon,klev,fluxt,"flux_vdf","Turbulent flux")
1353	call writeg1d(klon,klev,flux_ajs,"flux_ajs","Dry adjust. flux")
1354	call writeg1d(klon,klev,flux_ec,"flux_ec","Ec flux")
1355	c call writeg1d(klon,1,solsw,"surfsw","Net SW flux at surface")
1356	c call writeg1d(klon,1,sollw,"surflw","Net LW flux at surface")
1357	call writeg1d(klon,1,radsol,"surfnet","Net flux at surface")
1358	call writeg1d(klon,klev,d_t_vdf,"dt_vdf","DT from clmain")
1359	call writeg1d(klon,klev,d_t_ajs,"dt_ajs","DT from ajsec")
1360	call writeg1d(klon,klev,d_t_ec,"dt_ec","DT from Ec")
1361	endif
1362
1363	c====================================================================
1364	c Si c'est la fin, il faut conserver l'etat de redemarrage
1365	c====================================================================
1366	c
1367	IF (lafin) THEN
1368	itau_phy = itau_phy + itap
1369	c REMETTRE TOUS LES PARAMETRES POUR OROGW...
1370	CALL phyredem ("restartphy.nc",dtime,radpas,
1371	. rlatd, rlond, ftsol, ftsoil,
1372	. falbe,
1373	. solsw, sollw,dlw,
1374	. radsol,
1375	. zmea,zstd,zsig,zgam,zthe,zpic,zval,
1376	. t_ancien)
1377
1378	c--------------FLOTT
1379	CMODEB LOTT
1380	C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
1381	C DU BILAN DE MOMENT ANGULAIRE.
1382	if (bilansmc.eq.1) then
1383	write(,)'Fermeture de aaam_bud.out (FL Vous parle)'
1384	close(27)
1385	close(28)
1386	endif !bilansmc
1387	CMODFIN
1388	c-------------
1389	c--------------SLEBONNOIS
1390	C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
1391	C DES BALLONS
1392	if (ballons.eq.1) then
1393	write(,)'Fermeture des ballons*.out'
1394	close(30)
1395	close(31)
1396	close(32)
1397	close(33)
1398	close(34)
1399	endif !ballons
1400	c-------------
1401	ENDIF
1402
1403
1404	RETURN
1405	END
1406
1407
1408
1409	***********************************************************************
1410	***********************************************************************
1411	***********************************************************************
1412	***********************************************************************
1413	***********************************************************************
1414	***********************************************************************
1415	***********************************************************************
1416	***********************************************************************
1417
1418	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
1419	IMPLICIT none
1420	c
1421	c Tranformer une variable de la grille physique a
1422	c la grille d'ecriture
1423	c
1424	INTEGER nfield,nlon,iim,jjmp1, jjm
1425	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
1426	c
1427	INTEGER i, n, ig
1428	c
1429	jjm = jjmp1 - 1
1430	DO n = 1, nfield
1431	DO i=1,iim
1432	ecrit(i,n) = fi(1,n)
1433	ecrit(i+jjm*iim,n) = fi(nlon,n)
1434	ENDDO
1435	DO ig = 1, nlon - 2
1436	ecrit(iim+ig,n) = fi(1+ig,n)
1437	ENDDO
1438	ENDDO
1439	RETURN
1440	END

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

Download in other formats:

Original Format