Context Navigation

physiq.F @ 105

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

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

Download in other formats:

Original Format