Context Navigation

physiq.F @ 957

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

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

Download in other formats:

Original Format