Context Navigation

physiq.F @ 119

Last change on this file since 119 was 106, checked in by slebonnois, 14 years ago
SL : mise a jour de phytitan pour etre conforme aux sources actuelles utilisees sur gnome.
File size: 48.2 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 = 0
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	c zzlay(i,l)=zphi(i,l)/RG
751	c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
752	zzlay(i,l)=RGRARA/(RG*RA-zphi(i,l))-RA
753	ENDDO
754	ENDDO
755	DO i=1,klon
756	c zzlev(i,1)=0.
757	c CORRECTION 13/01/2011
758	c (correspond a la position de la surface en ce point vs RA)
759	zzlev(i,1)=pphis(i)/RG
760	ENDDO
761	DO l=2,klev
762	DO i=1,klon
763	z1=(pplay(i,l-1)+paprs(i,l))/(pplay(i,l-1)-paprs(i,l))
764	z2=(paprs(i,l)+pplay(i,l))/(paprs(i,l)-pplay(i,l))
765	zzlev(i,l)=(z1zzlay(i,l-1)+z2zzlay(i,l))/(z1+z2)
766	ENDDO
767	ENDDO
768	DO i=1,klon
769	zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
770	ENDDO
771
772	c- - - - - - - - - - - - - - - -
773	c DIAGNOSTIQUE GRILLE VERTICALE
774	c- - - - - - - - - - - - - - - -
775	c print*,"DIAGNOSTIQUE GRILLE VERTICALE"
776	c i=klon/2
777	c print*,"Niveau Pression Altitude (lev puis lay)"
778	c do l=1,klev
779	c print*,l,paprs(i,l),zzlev(i,l)
780	c print*,l,pplay(i,l),zzlay(i,l)
781	c enddo
782	c print*,klev+1,paprs(i,klev+1),zzlev(i,klev+1)
783	c stop
784
785	c====================================================================
786	c
787	c Verifier les temperatures
788	c
789	CALL hgardfou(t_seri,ftsol,'debutphy')
790	c====================================================================
791	c
792	c Incrementer le compteur de la physique
793	c
794	itap = itap + 1
795
796	c====================================================================
797	c
798	c Epaisseurs couches
799
800	DO k = 1, klev
801	DO i = 1, klon
802	delp(i,k) = paprs(i,k)-paprs(i,k+1)
803	ENDDO
804	ENDDO
805
806
807
808
809	c====================================================================
810	c ORBITE ET ECLAIREMENT
811	c====================================================================
812
813
814	c Pour TITAN:
815	c calcul de la longitude solaire
816	CALL solarlong(rjourvrai,zls)
817	print,'Ls',zls180./RPI ! zls est en radians !!
818	CALL orbite(zls,dist,pdecli)
819
820	c dans zenang, Ls en degres ; dans mucorr, Ls en radians
821	IF (cycle_diurne) THEN
822	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
823	CALL zenang(zls*180./RPI,gmtime,zdtime,rlatd,rlond,rmu0,fract)
824	ELSE
825	call mucorr(klon,zls,rlatd,rmu0,fract)
826	ENDIF
827
828
829
830
831	c====================================================================
832	c COUCHE LIMITE
833	c====================================================================
834
835	c-------------------------------
836	c TEST: on ne tient pas compte des calculs de clmain mais on force
837	c l'equilibre radiatif du sol
838	if (1.eq.0) then
839	if (debut) then
840	print*,"ATTENTION, CLMAIN SHUNTEE..."
841	endif
842
843	DO i = 1, klon
844	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
845	fder(i) = 0.0e0
846	dlw(i) = 0.0e0
847	ENDDO
848
849	c Incrementer la temperature du sol
850	c
851	DO i = 1, klon
852	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
853	ftsol(i) = ftsol(i) + d_ts(i)
854	do j=1,nsoilmx
855	ftsoil(i,j)=ftsol(i)
856	enddo
857	ENDDO
858
859	c-------------------------------
860	else
861	c-------------------------------
862
863	fder = dlw
864
865	c print*,"radsol avant clmain=",radsol(klon/2)
866	c print*,"solsw avant clmain=",solsw(klon/2)
867	c print*,"sollw avant clmain=",sollw(klon/2)
868
869	c CLMAIN
870
871	! ADAPTATION GCM POUR CP(T)
872	CALL clmain(dtime,itap,
873	e t_seri,u_seri,v_seri,
874	e rmu0,
875	e ftsol,
876	$ ftsoil,
877	$ paprs,pplay,ppk,radsol,falbe,
878	e solsw, sollw, sollwdown, fder,
879	e rlond, rlatd, cuphy, cvphy,
880	e debut, lafin,
881	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
882	s fluxt,fluxu,fluxv,cdragh,cdragm,
883	s dsens,
884	s ycoefh,yu1,yv1)
885
886	c print*,"radsol apres clmain=",radsol(klon/2)
887	c print*,"solsw apres clmain=",solsw(klon/2)
888	c print*,"sollw apres clmain=",sollw(klon/2)
889
890	CXXX Incrementation des flux
891	DO i = 1, klon
892	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
893	fder(i) = dlw(i) + dsens(i)
894	ENDDO
895	CXXX
896
897	DO k = 1, klev
898	DO i = 1, klon
899	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
900	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
901	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
902	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
903	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
904	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
905	ENDDO
906	ENDDO
907	c
908	c print,"d_t_vdf1=",d_t_vdf(1,:)dtime
909	c print,"d_t_vdf2=",d_t_vdf(klon/2,:)dtime
910	c print,"d_t_vdf3=",d_t_vdf(klon,:)dtime
911	c print,"d_u_vdf=",d_u_vdf(klon/2,:)dtime
912	c print,"d_v_vdf=",d_v_vdf(klon/2,:)dtime
913
914	C TRACEURS
915
916	d_tr_vdf = 0.
917	if (iflag_trac.eq.1) then
918	DO iq=1, nqmax
919	CALL cltrac(dtime,ycoefh,t_seri,
920	s tr_seri(1,1,iq), source,
921	e paprs, pplay, delp,
922	s d_tr_vdf(1,1,iq))
923
924	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
925	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
926	ENDDO
927	endif
928
929
930	IF (if_ebil.ge.2) THEN
931	ztit='after clmain'
932	CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
933	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
934	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
935	call diagphy(airephy,ztit,ip_ebil
936	e , zero_v, zero_v, zero_v, zero_v, sens
937	e , zero_v, zero_v, zero_v, ztsol
938	e , d_h_vcol, d_qt, d_ec
939	s , fs_bound, fq_bound )
940	END IF
941	C
942	c
943	c Incrementer la temperature du sol
944	c
945	c print*,'Tsol avant clmain:',ftsol(klon/2)
946	DO i = 1, klon
947	ftsol(i) = ftsol(i) + d_ts(i)
948	ENDDO
949	c print*,'DTsol apres clmain:',d_ts(klon/2)
950	c print*,'Tsol apres clmain:',ftsol(klon/2)
951
952	c Calculer la derive du flux infrarouge
953	c
954	DO i = 1, klon
955	dlw(i) = - 4.0RSIGMAftsol(i)**3
956	ENDDO
957
958	c-------------------------------
959	endif ! fin du TEST
960
961	c
962	c Appeler l'ajustement sec
963	c
964	c===================================================================
965	c Convection seche
966	c===================================================================
967	c
968	d_t_ajs(:,:)=0.
969	d_u_ajs(:,:)=0.
970	d_v_ajs(:,:)=0.
971	d_tr_ajs(:,:,:)=0.
972	c
973	IF(prt_level>9)WRITE(lunout,*)
974	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
975	s ,iflag_ajs
976
977	if(iflag_ajs.eq.0) then
978	c Rien
979	c ====
980	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
981
982	else if(iflag_ajs.eq.1) then
983
984	c Ajustement sec
985	c ==============
986	IF(prt_level>9)WRITE(lunout,*)'ajsec'
987
988	! ADAPTATION GCM POUR CP(T)
989	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
990	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
991
992	! ADAPTATION GCM POUR CP(T)
993	do i=1,klon
994	flux_ajs(i,1) = 0.0
995	do j=2,klev
996	flux_ajs(i,j) = flux_ajs(i,j-1)
997	. + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
998	. *delp(i,j-1)
999	enddo
1000	enddo
1001
1002	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
1003	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
1004	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
1005	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
1006	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
1007	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
1008	if (iflag_trac.eq.1) then
1009	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
1010	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
1011	endif
1012
1013	c print,"d_t_ajs1=",d_t_ajs(1,:)dtime
1014	c print,"d_t_ajs2=",d_t_ajs(klon/2,:)dtime
1015	c print,"d_t_ajs3=",d_t_ajs(klon,:)dtime
1016	c print,"d_u_ajs=",d_u_ajs(klon/2,:)dtime
1017	c print,"d_v_ajs=",d_v_ajs(klon/2,:)dtime
1018
1019	endif
1020	c
1021	IF (if_ebil.ge.2) THEN
1022	ztit='after dry_adjust'
1023	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1024	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1025	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1026	call diagphy(airephy,ztit,ip_ebil
1027	e , zero_v, zero_v, zero_v, zero_v, sens
1028	e , zero_v, zero_v, zero_v, ztsol
1029	e , d_h_vcol, d_qt, d_ec
1030	s , fs_bound, fq_bound )
1031	END IF
1032
1033	c====================================================================
1034	c MICROPHYSIQUE ET CHIMIE
1035	c====================================================================
1036
1037	d_tr_mph(:,:,:)=0.
1038	d_tr_kim(:,:,:)=0.
1039
1040	c on recupere tr_seri inchange, d_tr_micro, d_tr_chim, tous les trois sur nqmax
1041	c on recupere aussi qaer pour le mettre dans les sorties
1042	c si microfi=1, sortie de qaer(1:nmicro)
1043	c si nmicro != nqmax et si chimi, sortie de tr_seri(nmicro+1:nqmax)
1044
1045	c faire un test comme pour rayonnement, avec chimi en plus comme flag,
1046	c pour voir si chimie appelee -> bouleen, qui passe dans phytrac.
1047	c faut aussi le pas de temps chimique: dtimechim, a passer..
1048
1049	appel_chim = 0
1050	IF (MOD(itapchim,chimpas).EQ.0) THEN
1051	c print*,'CHIMIE ',
1052	c $ ' (itapchim=',itapchim,'/chimpas=',chimpas,')'
1053	appel_chim = 1
1054	itapchim = 0
1055	ENDIF
1056	itapchim = itapchim + 1
1057
1058
1059	if (iflag_trac.eq.1) then
1060	call phytrac (debut,lafin,
1061	. nqmax,nmicro,dtime,appel_chim,dtimechim,
1062	. paprs,pplay,delp,t,rmu0,fract,pdecli,zls,
1063	. tr_seri,qaer,d_tr_mph,d_tr_kim)
1064
1065	if (microfi.eq.1) then
1066	tr_seri(:,:,1:nmicro) = tr_seri(:,:,1:nmicro)
1067	. + d_tr_mph(:,:,1:nmicro)*dtime
1068	endif
1069	if ((chimi).and.(nqmax.gt.nmicro)) then
1070	c chimie:
1071	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_kim(:,:,:)*dtime
1072	c condensation:
1073	tr_seri(:,:,nmicro+1:nqmax) = tr_seri(:,:,nmicro+1:nqmax)
1074	. + d_tr_mph(:,:,nmicro+1:nqmax)*dtime
1075	endif
1076	endif
1077
1078	c------------------
1079	c test condensation
1080	c do i=1,nqmax
1081	c if(tname(i).eq."HCN") then
1082	c print*,"HCN="
1083	c do k=1,klev
1084	c print,k,tr_seri(klon/2,k,i),d_tr_mph(klon/2,k,i)dtime
1085	c v ,d_tr_kim(klon/2,k,i)*dtime
1086	c enddo
1087	c stop
1088	c endif
1089	c enddo
1090	c------------------
1091
1092	c====================================================================
1093	c RAYONNEMENT
1094	c====================================================================
1095
1096	IF (MOD(itaprad,radpas).EQ.0) THEN
1097	c print*,'RAYONNEMENT ',
1098	c $ ' (itaprad=',itaprad,'/radpas=',radpas,')'
1099
1100	c ATTENTION, (klon/2) ne marche pas toujours............
1101	c print*,"radsol avant radlwsw=",radsol(klon/2)
1102	c print*,"solsw avant radlwsw=",solsw(klon/2)
1103	c print*,"sollw avant radlwsw=",sollw(klon/2)
1104	c print*,"avant radlwsw"
1105	CALL radlwsw
1106	e (dist, rmu0, fract, dtimerad, zzlev,
1107	e paprs, pplay,ftsol, t_seri, nqmax, nmicro,
1108	c tr_seri, qaer,
1109	s heat,cool,radsol,
1110	s topsw,toplw,solsw,sollw,
1111	s sollwdown,
1112	s lwnet, swnet)
1113	c print*,"apres radlwsw"
1114
1115	c print*,"radsol apres radlwsw=",radsol(klon/2)
1116	c print*,"solsw apres radlwsw=",solsw(klon/2)
1117	c print*,"sollw apres radlwsw=",sollw(klon/2)
1118	itaprad = 0
1119	DO k = 1, klev
1120	DO i = 1, klon
1121	dtrad(i,k) = heat(i,k)-cool(i,k) !K/s
1122	ENDDO
1123	ENDDO
1124	c print*,"heat (K/s) =",heat(klon/2,:)
1125	c print*,"cool (K/s) =",cool(klon/2,:)
1126	c print*,"dtrad1 (K/s) =",dtrad(1,:)
1127	c print*,"dtrad2 (K/s) =",dtrad(klon/2,:)
1128	c print*,"dtrad3 (K/s) =",dtrad(klon,:)
1129
1130	ENDIF
1131	itaprad = itaprad + 1
1132	c====================================================================
1133	c
1134	c Ajouter la tendance des rayonnements (tous les pas)
1135	c
1136	DO k = 1, klev
1137	DO i = 1, klon
1138	t_seri(i,k) = t_seri(i,k) + dtrad(i,k) * dtime
1139	ENDDO
1140	ENDDO
1141
1142	IF (if_ebil.ge.2) THEN
1143	ztit='after rad'
1144	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1145	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1146	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1147	call diagphy(airephy,ztit,ip_ebil
1148	e , topsw, toplw, solsw, sollw, zero_v
1149	e , zero_v, zero_v, zero_v, ztsol
1150	e , d_h_vcol, d_qt, d_ec
1151	s , fs_bound, fq_bound )
1152	END IF
1153	c
1154
1155	c====================================================================
1156	c+jld ec_conser
1157	DO k = 1, klev
1158	DO i = 1, klon
1159	d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
1160	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
1161	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1162	d_t_ec(i,k) = d_t_ec(i,k)/dtime
1163	END DO
1164	END DO
1165	do i=1,klon
1166	flux_ec(i,1) = 0.0
1167	do j=2,klev
1168	flux_ec(i,j) = flux_ec(i,j-1)
1169	. +cpdet(t_seri(i,j-1))/RGd_t_ec(i,j-1)delp(i,j-1)
1170	enddo
1171	enddo
1172
1173	c-jld ec_conser
1174	c====================================================================
1175
1176	IF (if_ebil.ge.1) THEN
1177	ztit='after physic'
1178	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
1179	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1180	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1181	C Comme les tendances de la physique sont ajoute dans la dynamique,
1182	C on devrait avoir que la variation d'entalpie par la dynamique
1183	C est egale a la variation de la physique au pas de temps precedent.
1184	C Donc la somme de ces 2 variations devrait etre nulle.
1185	call diagphy(airephy,ztit,ip_ebil
1186	e , topsw, toplw, solsw, sollw, sens
1187	e , zero_v, zero_v, zero_v, ztsol
1188	e , d_h_vcol, d_qt, d_ec
1189	s , fs_bound, fq_bound )
1190	C
1191	d_h_vcol_phy=d_h_vcol
1192	C
1193	END IF
1194	C
1195	c====================================================================
1196	c Calcul des gravity waves FLOTT
1197	c====================================================================
1198	c
1199	c if (ok_orodr.or.ok_gw_nonoro) then
1200	cc CALCUL DE N2
1201	c do i=1,klon
1202	c do k=2,klev
1203	c ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1204	c zpklev(i,k) = sqrt(ppk(i,k)*ppk(i,k-1))
1205	c enddo
1206	c enddo
1207	c call t2tpot(klon*klev,ztlev, ztetalev,zpklev)
1208	c call t2tpot(klon*klev,t_seri,ztetalay,ppk)
1209	c do i=1,klon
1210	c do k=2,klev
1211	c zdtetalev(i,k) = ztetalay(i,k)-ztetalay(i,k-1)
1212	c zdzlev(i,k) = (zphi(i,k)-zphi(i,k-1))/RG
1213	c zn2(i,k) = RGzdtetalev(i,k)/(ztetalev(i,k)zdzlev(i,k))
1214	c zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
1215	c enddo
1216	c enddo
1217	c
1218	c endif
1219	c
1220	cc ----------------------------ORODRAG
1221	c IF (ok_orodr) THEN
1222	cc
1223	cc selection des points pour lesquels le shema est actif:
1224	c igwd=0
1225	c DO i=1,klon
1226	c itest(i)=0
1227	cc IF ((zstd(i).gt.10.0)) THEN
1228	c IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1229	c itest(i)=1
1230	c igwd=igwd+1
1231	c idx(igwd)=i
1232	c ENDIF
1233	c ENDDO
1234	cc igwdim=MAX(1,igwd)
1235	cc
1236	cc A ADAPTER POUR TITAN !!!
1237	c CALL drag_noro(klon,klev,dtime,paprs,pplay,zphi,zn2,
1238	c e zmea,zstd, zsig, zgam, zthe,zpic,zval,
1239	c e igwd,idx,itest,
1240	c e t_seri, u_seri, v_seri,
1241	c s zulow, zvlow, zustrdr, zvstrdr,
1242	c s d_t_oro, d_u_oro, d_v_oro)
1243	c
1244	cc ajout des tendances
1245	c t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
1246	c d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
1247	c u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
1248	c d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
1249	c v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
1250	c d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
1251	cc
1252	c ELSE
1253	c d_t_oro = 0.
1254	c d_u_oro = 0.
1255	c d_v_oro = 0.
1256	c zustrdr = 0.
1257	c zvstrdr = 0.
1258	cc
1259	c ENDIF ! fin de test sur ok_orodr
1260	cc
1261	cc ----------------------------OROLIFT
1262	c IF (ok_orolf) THEN
1263	cc
1264	cc selection des points pour lesquels le shema est actif:
1265	c igwd=0
1266	c DO i=1,klon
1267	c itest(i)=0
1268	c IF ((zpic(i)-zmea(i)).GT.100.) THEN
1269	c itest(i)=1
1270	c igwd=igwd+1
1271	c idx(igwd)=i
1272	c ENDIF
1273	c ENDDO
1274	cc igwdim=MAX(1,igwd)
1275	cc
1276	cc A ADAPTER POUR VENUS!!!
1277	cc CALL lift_noro(klon,klev,dtime,paprs,pplay,
1278	cc e rlatd,zmea,zstd,zpic,zgam,zthe,zpic,zval,
1279	cc e igwd,idx,itest,
1280	cc e t_seri, u_seri, v_seri,
1281	cc s zulow, zvlow, zustrli, zvstrli,
1282	cc s d_t_lif, d_u_lif, d_v_lif )
1283	c
1284	cc
1285	cc ajout des tendances
1286	c t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
1287	c d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
1288	c u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
1289	c d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
1290	c v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
1291	c d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
1292	cc
1293	c ELSE
1294	c d_t_lif = 0.
1295	c d_u_lif = 0.
1296	c d_v_lif = 0.
1297	c zustrli = 0.
1298	c zvstrli = 0.
1299	cc
1300	c ENDIF ! fin de test sur ok_orolf
1301	c
1302	cc ---------------------------- NON-ORO GRAVITY WAVES
1303	c IF(ok_gw_nonoro) then
1304	c
1305	c call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
1306	c e t_seri, u_seri, v_seri,
1307	c o zustrhi,zvstrhi,
1308	c o d_t_hin, d_u_hin, d_v_hin)
1309	c
1310	cc ajout des tendances
1311	c
1312	c t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
1313	c d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
1314	c u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
1315	c d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
1316	c v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
1317	c d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
1318	c
1319	c ELSE
1320	c d_t_hin = 0.
1321	c d_u_hin = 0.
1322	c d_v_hin = 0.
1323	c zustrhi = 0.
1324	c zvstrhi = 0.
1325	c
1326	c ENDIF ! fin de test sur ok_gw_nonoro
1327	c
1328	c====================================================================
1329	c Transport de ballons
1330	c====================================================================
1331	if (ballons.eq.1) then
1332	CALL ballon(30,pdtphys,rjourvrai,gmtime,rlatd,rlond,
1333	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
1334	C t,pplay,u,v,zphi) ! alt above planet average radius
1335	endif !ballons
1336
1337	c====================================================================
1338	c Bilan de mmt angulaire
1339	c====================================================================
1340	if (bilansmc.eq.1) then
1341	CMODDEB FLOTT
1342	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
1343	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
1344
1345	DO i = 1, klon
1346	zustrph(i)=0.
1347	zvstrph(i)=0.
1348	zustrcl(i)=0.
1349	zvstrcl(i)=0.
1350	ENDDO
1351	DO k = 1, klev
1352	DO i = 1, klon
1353	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
1354	c (paprs(i,k)-paprs(i,k+1))/rg
1355	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
1356	c (paprs(i,k)-paprs(i,k+1))/rg
1357	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
1358	c (paprs(i,k)-paprs(i,k+1))/rg
1359	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
1360	c (paprs(i,k)-paprs(i,k+1))/rg
1361	ENDDO
1362	ENDDO
1363
1364	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
1365	C ra,rg,romega,
1366	C rlatd,rlond,pphis,
1367	C zustrdr,zustrli,zustrcl,
1368	C zvstrdr,zvstrli,zvstrcl,
1369	C paprs,u,v)
1370
1371	CCMODFIN FLOTT
1372	endif !bilansmc
1373
1374	c=======================================================================
1375	c SORTIES
1376	c=======================================================================
1377
1378	c Convertir les incrementations en tendances
1379	c
1380	DO k = 1, klev
1381	DO i = 1, klon
1382	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
1383	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
1384	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
1385	ENDDO
1386	ENDDO
1387	c print*,"vnatphy=",v(705,:)
1388	c print*,"unatphy=",u(705,:)
1389	c
1390	DO iq = 1, nqmax
1391	DO k = 1, klev
1392	DO i = 1, klon
1393	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
1394	ENDDO
1395	ENDDO
1396	ENDDO
1397
1398	c------------------------
1399	c Calcul moment cinetique
1400	c------------------------
1401	c TEST
1402	c mangtot = 0.0
1403	c DO k = 1, klev
1404	c DO i = 1, klon
1405	c mang(i,k) = RAcos(rlatd(i)RPI/180.)
1406	c . (u_seri(i,k)+RAcos(rlatd(i)RPI/180.)ROMEGA)
1407	c . airephy(i)delp(i,k)/RG
1408	c mangtot=mangtot+mang(i,k)
1409	c ENDDO
1410	c ENDDO
1411	c print*,"Moment cinetique total = ",mangtot
1412	c
1413	c------------------------
1414	c
1415	c Sauvegarder les valeurs de t et u a la fin de la physique:
1416	c
1417	DO k = 1, klev
1418	DO i = 1, klon
1419	u_ancien(i,k) = u_seri(i,k)
1420	t_ancien(i,k) = t_seri(i,k)
1421	ENDDO
1422	ENDDO
1423	c
1424	c incrementation de la tendance (repassee en extensif) sur qaer pour sorties
1425	if (microfi.eq.1) then
1426	do iq=1,nmicro
1427	DO l=1,llm
1428	DO i = 1, klon
1429	qaer(i,l,iq) = qaer(i,l,iq) +
1430	. (d_tr_mph(i,l,iq)delp(i,l)/RG)dtime
1431	ENDDO
1432	ENDDO
1433	enddo
1434	endif ! microfi
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