Context Navigation

physiq.F @ 91

Last change on this file since 91 was 3, checked in by slebonnois, 15 years ago

Creation de repertoires:

chantiers : pour communiquer sur nos projets de modifs
documentation : pour stocker les docs

Ajout de:

libf/phytitan : physique de Titan
libf/chimtitan: chimie de Titan
libf/phyvenus : physique de Venus

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

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

Download in other formats:

Original Format