Context Navigation

physiq.F @ 1017

Last change on this file since 1017 was 1017, checked in by emillour, 11 years ago

Common dynamics: (and collateral adaptations in Venus physics)
Improved cpdet routines in and additional sponge mode:

Additionnal sponge mode (trigered with "callsponge" flag), in line with the one used in the Generic and Martian GCM. This sponge is called whenever there is a dissipation step.
Improvement of the cpdet routines : created routines tpot2t_glo_p and t2tpot_glo_p which handle fields on the whole dynamics (scaler) grid, which are more efficient than calling tpot2t_p or t2tpot_p with slabs of data (generated use of intermediate copies of these chunks of data at every call)
Turned cpdet.F into a module cpdet_mod.F90 (and correspondingly adapted all routines in the Venus physics).

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

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

Download in other formats:

Original Format