Context Navigation

physiq.F @ 124

Last change on this file since 124 was 123, checked in by slebonnois, 14 years ago
SL: encore un bug oublie dans physiq.F
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-NBjours
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+gmtime,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	c====================================================================
829	c COUCHE LIMITE
830	c====================================================================
831
832	c-------------------------------
833	c TEST: on ne tient pas compte des calculs de clmain mais on force
834	c l'equilibre radiatif du sol
835	if (1.eq.0) then
836	if (debut) then
837	print*,"ATTENTION, CLMAIN SHUNTEE..."
838	endif
839
840	DO i = 1, klon
841	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
842	fder(i) = 0.0e0
843	dlw(i) = 0.0e0
844	ENDDO
845
846	c Incrementer la temperature du sol
847	c
848	DO i = 1, klon
849	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
850	ftsol(i) = ftsol(i) + d_ts(i)
851	do j=1,nsoilmx
852	ftsoil(i,j)=ftsol(i)
853	enddo
854	ENDDO
855
856	c-------------------------------
857	else
858	c-------------------------------
859
860	fder = dlw
861
862	c print*,"radsol avant clmain=",radsol(klon/2)
863	c print*,"solsw avant clmain=",solsw(klon/2)
864	c print*,"sollw avant clmain=",sollw(klon/2)
865
866	c CLMAIN
867
868	! ADAPTATION GCM POUR CP(T)
869	CALL clmain(dtime,itap,
870	e t_seri,u_seri,v_seri,
871	e rmu0,
872	e ftsol,
873	$ ftsoil,
874	$ paprs,pplay,ppk,radsol,falbe,
875	e solsw, sollw, sollwdown, fder,
876	e rlond, rlatd, cuphy, cvphy,
877	e debut, lafin,
878	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
879	s fluxt,fluxu,fluxv,cdragh,cdragm,
880	s dsens,
881	s ycoefh,yu1,yv1)
882
883	c print*,"radsol apres clmain=",radsol(klon/2)
884	c print*,"solsw apres clmain=",solsw(klon/2)
885	c print*,"sollw apres clmain=",sollw(klon/2)
886
887	CXXX Incrementation des flux
888	DO i = 1, klon
889	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
890	fder(i) = dlw(i) + dsens(i)
891	ENDDO
892	CXXX
893
894	DO k = 1, klev
895	DO i = 1, klon
896	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
897	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
898	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
899	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
900	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
901	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
902	ENDDO
903	ENDDO
904
905	c print,"d_t_vdf1=",d_t_vdf(1,:)dtime
906	c print,"d_t_vdf2=",d_t_vdf(klon/2,:)dtime
907	c print,"d_t_vdf3=",d_t_vdf(klon,:)dtime
908	c print,"d_u_vdf=",d_u_vdf(klon/2,:)dtime
909	c print,"d_v_vdf=",d_v_vdf(klon/2,:)dtime
910
911	C TRACEURS
912
913	d_tr_vdf = 0.
914	if (iflag_trac.eq.1) then
915	DO iq=1, nqmax
916	CALL cltrac(dtime,ycoefh,t_seri,
917	s tr_seri(1,1,iq), source,
918	e paprs, pplay, delp,
919	s d_tr_vdf(1,1,iq))
920
921	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
922	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
923	ENDDO
924	endif
925
926	IF (if_ebil.ge.2) THEN
927	ztit='after clmain'
928	CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
929	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
930	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
931	call diagphy(airephy,ztit,ip_ebil
932	e , zero_v, zero_v, zero_v, zero_v, sens
933	e , zero_v, zero_v, zero_v, ztsol
934	e , d_h_vcol, d_qt, d_ec
935	s , fs_bound, fq_bound )
936	END IF
937	C
938	c
939	c Incrementer la temperature du sol
940	c
941	c print*,'Tsol avant clmain:',ftsol(klon/2)
942	DO i = 1, klon
943	ftsol(i) = ftsol(i) + d_ts(i)
944	ENDDO
945	c print*,'DTsol apres clmain:',d_ts(klon/2)
946	c print*,'Tsol apres clmain:',ftsol(klon/2)
947
948	c Calculer la derive du flux infrarouge
949	c
950	DO i = 1, klon
951	dlw(i) = - 4.0RSIGMAftsol(i)**3
952	ENDDO
953
954	c-------------------------------
955	endif ! fin du TEST
956
957	c
958	c Appeler l'ajustement sec
959	c
960	c===================================================================
961	c Convection seche
962	c===================================================================
963	c
964	d_t_ajs(:,:)=0.
965	d_u_ajs(:,:)=0.
966	d_v_ajs(:,:)=0.
967	d_tr_ajs(:,:,:)=0.
968	c
969	IF(prt_level>9)WRITE(lunout,*)
970	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
971	s ,iflag_ajs
972
973	if(iflag_ajs.eq.0) then
974	c Rien
975	c ====
976	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
977
978	else if(iflag_ajs.eq.1) then
979
980	c Ajustement sec
981	c ==============
982	IF(prt_level>9)WRITE(lunout,*)'ajsec'
983
984	! ADAPTATION GCM POUR CP(T)
985	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
986	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
987
988	! ADAPTATION GCM POUR CP(T)
989	do i=1,klon
990	flux_ajs(i,1) = 0.0
991	do j=2,klev
992	flux_ajs(i,j) = flux_ajs(i,j-1)
993	. + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
994	. *delp(i,j-1)
995	enddo
996	enddo
997
998	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
999	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
1000	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
1001	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
1002	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
1003	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
1004	if (iflag_trac.eq.1) then
1005	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
1006	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
1007	endif
1008
1009	c print,"d_t_ajs1=",d_t_ajs(1,:)dtime
1010	c print,"d_t_ajs2=",d_t_ajs(klon/2,:)dtime
1011	c print,"d_t_ajs3=",d_t_ajs(klon,:)dtime
1012	c print,"d_u_ajs=",d_u_ajs(klon/2,:)dtime
1013	c print,"d_v_ajs=",d_v_ajs(klon/2,:)dtime
1014
1015	endif
1016	c
1017	IF (if_ebil.ge.2) THEN
1018	ztit='after dry_adjust'
1019	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1020	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1021	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1022	call diagphy(airephy,ztit,ip_ebil
1023	e , zero_v, zero_v, zero_v, zero_v, sens
1024	e , zero_v, zero_v, zero_v, ztsol
1025	e , d_h_vcol, d_qt, d_ec
1026	s , fs_bound, fq_bound )
1027	END IF
1028
1029	c====================================================================
1030	c MICROPHYSIQUE ET CHIMIE
1031	c====================================================================
1032
1033	d_tr_mph(:,:,:)=0.
1034	d_tr_kim(:,:,:)=0.
1035
1036	c on recupere tr_seri inchange, d_tr_micro, d_tr_chim, tous les trois sur nqmax
1037	c on recupere aussi qaer pour le mettre dans les sorties
1038	c si microfi=1, sortie de qaer(1:nmicro)
1039	c si nmicro != nqmax et si chimi, sortie de tr_seri(nmicro+1:nqmax)
1040
1041	c faire un test comme pour rayonnement, avec chimi en plus comme flag,
1042	c pour voir si chimie appelee -> bouleen, qui passe dans phytrac.
1043	c faut aussi le pas de temps chimique: dtimechim, a passer..
1044
1045	appel_chim = 0
1046	IF (MOD(itapchim,chimpas).EQ.0) THEN
1047	c print*,'CHIMIE ',
1048	c $ ' (itapchim=',itapchim,'/chimpas=',chimpas,')'
1049	appel_chim = 1
1050	itapchim = 0
1051	ENDIF
1052	itapchim = itapchim + 1
1053
1054	if (iflag_trac.eq.1) then
1055	call phytrac (debut,lafin,
1056	. nqmax,nmicro,dtime,appel_chim,dtimechim,
1057	. paprs,pplay,delp,t,rmu0,fract,pdecli,zls,
1058	. tr_seri,qaer,d_tr_mph,d_tr_kim)
1059
1060	if (microfi.eq.1) then
1061	tr_seri(:,:,1:nmicro) = tr_seri(:,:,1:nmicro)
1062	. + d_tr_mph(:,:,1:nmicro)*dtime
1063	endif
1064	if ((chimi).and.(nqmax.gt.nmicro)) then
1065	c chimie:
1066	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_kim(:,:,:)*dtime
1067	c condensation:
1068	tr_seri(:,:,nmicro+1:nqmax) = tr_seri(:,:,nmicro+1:nqmax)
1069	. + d_tr_mph(:,:,nmicro+1:nqmax)*dtime
1070	endif
1071
1072	endif
1073
1074	c------------------
1075	c test condensation
1076	c do i=1,nqmax
1077	c if(tname(i).eq."HCN") then
1078	c print*,"HCN="
1079	c do k=1,klev
1080	c print,k,tr_seri(klon/2,k,i),d_tr_mph(klon/2,k,i)dtime
1081	c v ,d_tr_kim(klon/2,k,i)*dtime
1082	c enddo
1083	c stop
1084	c endif
1085	c enddo
1086	c------------------
1087
1088	c====================================================================
1089	c RAYONNEMENT
1090	c====================================================================
1091
1092	IF (MOD(itaprad,radpas).EQ.0) THEN
1093	c print*,'RAYONNEMENT ',
1094	c $ ' (itaprad=',itaprad,'/radpas=',radpas,')'
1095
1096	c ATTENTION, (klon/2) ne marche pas toujours............
1097	c print*,"radsol avant radlwsw=",radsol(klon/2)
1098	c print*,"solsw avant radlwsw=",solsw(klon/2)
1099	c print*,"sollw avant radlwsw=",sollw(klon/2)
1100	c print*,"avant radlwsw"
1101	CALL radlwsw
1102	e (dist, rmu0, fract, dtimerad, zzlev,
1103	e paprs, pplay,ftsol, t_seri, nqmax, nmicro,
1104	c tr_seri, qaer,
1105	s heat,cool,radsol,
1106	s topsw,toplw,solsw,sollw,
1107	s sollwdown,
1108	s lwnet, swnet)
1109	c print*,"apres radlwsw"
1110
1111	c print*,"radsol apres radlwsw=",radsol(klon/2)
1112	c print*,"solsw apres radlwsw=",solsw(klon/2)
1113	c print*,"sollw apres radlwsw=",sollw(klon/2)
1114	itaprad = 0
1115	DO k = 1, klev
1116	DO i = 1, klon
1117	dtrad(i,k) = heat(i,k)-cool(i,k) !K/s
1118	ENDDO
1119	ENDDO
1120	c print*,"heat (K/s) =",heat(klon/2,:)
1121	c print*,"cool (K/s) =",cool(klon/2,:)
1122	c print*,"dtrad1 (K/s) =",dtrad(1,:)
1123	c print*,"dtrad2 (K/s) =",dtrad(klon/2,:)
1124	c print*,"dtrad3 (K/s) =",dtrad(klon,:)
1125
1126	ENDIF
1127	itaprad = itaprad + 1
1128	c====================================================================
1129	c
1130	c Ajouter la tendance des rayonnements (tous les pas)
1131	c
1132	DO k = 1, klev
1133	DO i = 1, klon
1134	t_seri(i,k) = t_seri(i,k) + dtrad(i,k) * dtime
1135	ENDDO
1136	ENDDO
1137
1138	IF (if_ebil.ge.2) THEN
1139	ztit='after rad'
1140	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1141	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1142	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1143	call diagphy(airephy,ztit,ip_ebil
1144	e , topsw, toplw, solsw, sollw, zero_v
1145	e , zero_v, zero_v, zero_v, ztsol
1146	e , d_h_vcol, d_qt, d_ec
1147	s , fs_bound, fq_bound )
1148	END IF
1149	c
1150
1151	c====================================================================
1152	c+jld ec_conser
1153	DO k = 1, klev
1154	DO i = 1, klon
1155	d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
1156	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
1157	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1158	d_t_ec(i,k) = d_t_ec(i,k)/dtime
1159	END DO
1160	END DO
1161	do i=1,klon
1162	flux_ec(i,1) = 0.0
1163	do j=2,klev
1164	flux_ec(i,j) = flux_ec(i,j-1)
1165	. +cpdet(t_seri(i,j-1))/RGd_t_ec(i,j-1)delp(i,j-1)
1166	enddo
1167	enddo
1168
1169	c-jld ec_conser
1170	c====================================================================
1171
1172	IF (if_ebil.ge.1) THEN
1173	ztit='after physic'
1174	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
1175	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1176	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1177	C Comme les tendances de la physique sont ajoute dans la dynamique,
1178	C on devrait avoir que la variation d'entalpie par la dynamique
1179	C est egale a la variation de la physique au pas de temps precedent.
1180	C Donc la somme de ces 2 variations devrait etre nulle.
1181	call diagphy(airephy,ztit,ip_ebil
1182	e , topsw, toplw, solsw, sollw, sens
1183	e , zero_v, zero_v, zero_v, ztsol
1184	e , d_h_vcol, d_qt, d_ec
1185	s , fs_bound, fq_bound )
1186	C
1187	d_h_vcol_phy=d_h_vcol
1188	C
1189	END IF
1190	C
1191	c====================================================================
1192	c Calcul des gravity waves FLOTT
1193	c====================================================================
1194	c
1195	c if (ok_orodr.or.ok_gw_nonoro) then
1196	cc CALCUL DE N2
1197	c do i=1,klon
1198	c do k=2,klev
1199	c ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1200	c zpklev(i,k) = sqrt(ppk(i,k)*ppk(i,k-1))
1201	c enddo
1202	c enddo
1203	c call t2tpot(klon*klev,ztlev, ztetalev,zpklev)
1204	c call t2tpot(klon*klev,t_seri,ztetalay,ppk)
1205	c do i=1,klon
1206	c do k=2,klev
1207	c zdtetalev(i,k) = ztetalay(i,k)-ztetalay(i,k-1)
1208	c zdzlev(i,k) = (zphi(i,k)-zphi(i,k-1))/RG
1209	c zn2(i,k) = RGzdtetalev(i,k)/(ztetalev(i,k)zdzlev(i,k))
1210	c zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
1211	c enddo
1212	c enddo
1213	c
1214	c endif
1215	c
1216	cc ----------------------------ORODRAG
1217	c IF (ok_orodr) THEN
1218	cc
1219	cc selection des points pour lesquels le shema est actif:
1220	c igwd=0
1221	c DO i=1,klon
1222	c itest(i)=0
1223	cc IF ((zstd(i).gt.10.0)) THEN
1224	c IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1225	c itest(i)=1
1226	c igwd=igwd+1
1227	c idx(igwd)=i
1228	c ENDIF
1229	c ENDDO
1230	cc igwdim=MAX(1,igwd)
1231	cc
1232	cc A ADAPTER POUR TITAN !!!
1233	c CALL drag_noro(klon,klev,dtime,paprs,pplay,zphi,zn2,
1234	c e zmea,zstd, zsig, zgam, zthe,zpic,zval,
1235	c e igwd,idx,itest,
1236	c e t_seri, u_seri, v_seri,
1237	c s zulow, zvlow, zustrdr, zvstrdr,
1238	c s d_t_oro, d_u_oro, d_v_oro)
1239	c
1240	cc ajout des tendances
1241	c t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
1242	c d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
1243	c u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
1244	c d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
1245	c v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
1246	c d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
1247	cc
1248	c ELSE
1249	c d_t_oro = 0.
1250	c d_u_oro = 0.
1251	c d_v_oro = 0.
1252	c zustrdr = 0.
1253	c zvstrdr = 0.
1254	cc
1255	c ENDIF ! fin de test sur ok_orodr
1256	cc
1257	cc ----------------------------OROLIFT
1258	c IF (ok_orolf) THEN
1259	cc
1260	cc selection des points pour lesquels le shema est actif:
1261	c igwd=0
1262	c DO i=1,klon
1263	c itest(i)=0
1264	c IF ((zpic(i)-zmea(i)).GT.100.) THEN
1265	c itest(i)=1
1266	c igwd=igwd+1
1267	c idx(igwd)=i
1268	c ENDIF
1269	c ENDDO
1270	cc igwdim=MAX(1,igwd)
1271	cc
1272	cc A ADAPTER POUR VENUS!!!
1273	cc CALL lift_noro(klon,klev,dtime,paprs,pplay,
1274	cc e rlatd,zmea,zstd,zpic,zgam,zthe,zpic,zval,
1275	cc e igwd,idx,itest,
1276	cc e t_seri, u_seri, v_seri,
1277	cc s zulow, zvlow, zustrli, zvstrli,
1278	cc s d_t_lif, d_u_lif, d_v_lif )
1279	c
1280	cc
1281	cc ajout des tendances
1282	c t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
1283	c d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
1284	c u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
1285	c d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
1286	c v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
1287	c d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
1288	cc
1289	c ELSE
1290	c d_t_lif = 0.
1291	c d_u_lif = 0.
1292	c d_v_lif = 0.
1293	c zustrli = 0.
1294	c zvstrli = 0.
1295	cc
1296	c ENDIF ! fin de test sur ok_orolf
1297	c
1298	cc ---------------------------- NON-ORO GRAVITY WAVES
1299	c IF(ok_gw_nonoro) then
1300	c
1301	c call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
1302	c e t_seri, u_seri, v_seri,
1303	c o zustrhi,zvstrhi,
1304	c o d_t_hin, d_u_hin, d_v_hin)
1305	c
1306	cc ajout des tendances
1307	c
1308	c t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
1309	c d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
1310	c u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
1311	c d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
1312	c v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
1313	c d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
1314	c
1315	c ELSE
1316	c d_t_hin = 0.
1317	c d_u_hin = 0.
1318	c d_v_hin = 0.
1319	c zustrhi = 0.
1320	c zvstrhi = 0.
1321	c
1322	c ENDIF ! fin de test sur ok_gw_nonoro
1323	c
1324	c====================================================================
1325	c Transport de ballons
1326	c====================================================================
1327	if (ballons.eq.1) then
1328	CALL ballon(30,pdtphys,rjourvrai,gmtime,rlatd,rlond,
1329	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
1330	C t,pplay,u,v,zphi) ! alt above planet average radius
1331	endif !ballons
1332
1333	c====================================================================
1334	c Bilan de mmt angulaire
1335	c====================================================================
1336	if (bilansmc.eq.1) then
1337	CMODDEB FLOTT
1338	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
1339	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
1340
1341	DO i = 1, klon
1342	zustrph(i)=0.
1343	zvstrph(i)=0.
1344	zustrcl(i)=0.
1345	zvstrcl(i)=0.
1346	ENDDO
1347	DO k = 1, klev
1348	DO i = 1, klon
1349	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
1350	c (paprs(i,k)-paprs(i,k+1))/rg
1351	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
1352	c (paprs(i,k)-paprs(i,k+1))/rg
1353	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
1354	c (paprs(i,k)-paprs(i,k+1))/rg
1355	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
1356	c (paprs(i,k)-paprs(i,k+1))/rg
1357	ENDDO
1358	ENDDO
1359
1360	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
1361	C ra,rg,romega,
1362	C rlatd,rlond,pphis,
1363	C zustrdr,zustrli,zustrcl,
1364	C zvstrdr,zvstrli,zvstrcl,
1365	C paprs,u,v)
1366
1367	CCMODFIN FLOTT
1368	endif !bilansmc
1369
1370	c=======================================================================
1371	c SORTIES
1372	c=======================================================================
1373
1374	c Convertir les incrementations en tendances
1375	c
1376	DO k = 1, klev
1377	DO i = 1, klon
1378	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
1379	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
1380	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
1381	ENDDO
1382	ENDDO
1383	c print*,"vnatphy=",v(705,:)
1384	c print*,"unatphy=",u(705,:)
1385	c
1386	DO iq = 1, nqmax
1387	DO k = 1, klev
1388	DO i = 1, klon
1389	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
1390	ENDDO
1391	ENDDO
1392	ENDDO
1393
1394	c------------------------
1395	c Calcul moment cinetique
1396	c------------------------
1397	c TEST
1398	c mangtot = 0.0
1399	c DO k = 1, klev
1400	c DO i = 1, klon
1401	c mang(i,k) = RAcos(rlatd(i)RPI/180.)
1402	c . (u_seri(i,k)+RAcos(rlatd(i)RPI/180.)ROMEGA)
1403	c . airephy(i)delp(i,k)/RG
1404	c mangtot=mangtot+mang(i,k)
1405	c ENDDO
1406	c ENDDO
1407	c print*,"Moment cinetique total = ",mangtot
1408	c
1409	c------------------------
1410	c
1411	c Sauvegarder les valeurs de t et u a la fin de la physique:
1412	c
1413	DO k = 1, klev
1414	DO i = 1, klon
1415	u_ancien(i,k) = u_seri(i,k)
1416	t_ancien(i,k) = t_seri(i,k)
1417	ENDDO
1418	ENDDO
1419	c
1420	c incrementation de la tendance (repassee en extensif) sur qaer pour sorties
1421	if (microfi.eq.1) then
1422	do iq=1,nmicro
1423	DO l=1,llm
1424	DO i = 1, klon
1425	qaer(i,l,iq) = qaer(i,l,iq) +
1426	. (d_tr_mph(i,l,iq)delp(i,l)/RG)dtime
1427	ENDDO
1428	ENDDO
1429	enddo
1430	endif ! microfi
1431	c=============================================================
1432	c Ecriture des sorties
1433	c=============================================================
1434
1435	#ifdef CPP_IOIPSL
1436
1437	#ifdef histmth
1438	#include "write_histmth.h"
1439	#endif
1440
1441	#ifdef histday
1442	#include "write_histday.h"
1443	#endif
1444
1445	#ifdef histins
1446	#include "write_histins.h"
1447	#endif
1448
1449	#endif
1450
1451	c ---------- Sorties testphys1d -------------
1452
1453	if (klon.eq.1) then
1454	call writeg1d(klon,klev,t_seri,"Temp","Temperature")
1455	call writeg1d(klon,1,ftsol,"tsurf","Surface Temp")
1456	DO k = 1, klev
1457	DO i = 1, klon
1458	c tmpout(i,k) = real(heat(i,k))
1459	tmpout(i,k) = heat(i,k)
1460	ENDDO
1461	ENDDO
1462	call writeg1d(klon,klev,tmpout,
1463	. "heat","Solar heating")
1464	DO k = 1, klev
1465	DO i = 1, klon
1466	c tmpout(i,k) = real(dtrad(i,k))
1467	tmpout(i,k) = dtrad(i,k)
1468	ENDDO
1469	ENDDO
1470	call writeg1d(klon,klev,tmpout,
1471	. "dtrad","IR cooling")
1472	call writeg1d(klon,klev,lwnet,"lwnet","Net LW flux")
1473	call writeg1d(klon,klev,swnet,"swnet","Net SW flux")
1474	call writeg1d(klon,klev,fluxt,"flux_vdf","Turbulent flux")
1475	call writeg1d(klon,klev,flux_ajs,"flux_ajs","Dry adjust. flux")
1476	call writeg1d(klon,klev,flux_ec,"flux_ec","Ec flux")
1477	c call writeg1d(klon,1,solsw,"surfsw","Net SW flux at surface")
1478	c call writeg1d(klon,1,sollw,"surflw","Net LW flux at surface")
1479	call writeg1d(klon,1,radsol,"surfnet","Net flux at surface")
1480	call writeg1d(klon,klev,d_t_vdf,"dt_vdf","DT from clmain")
1481	call writeg1d(klon,klev,d_t_ajs,"dt_ajs","DT from ajsec")
1482	call writeg1d(klon,klev,d_t_ec,"dt_ec","DT from Ec")
1483	endif
1484
1485	c====================================================================
1486	c Si c'est la fin, il faut conserver l'etat de redemarrage
1487	c====================================================================
1488	c
1489	IF (lafin) THEN
1490	itau_phy = itau_phy + itap
1491	c REMETTRE TOUS LES PARAMETRES POUR OROGW... A FAIRE POUR TITAN
1492	CALL phyredem ("restartphy.nc",dtime,radpas,chimpas,
1493	. rlatd, rlond, ftsol, ftsoil,
1494	. falbe,
1495	. solsw, sollw,dlw,
1496	. radsol,
1497	c . zmea,zstd,zsig,zgam,zthe,zpic,zval,
1498	. t_ancien)
1499
1500	c--------------FLOTT
1501	CMODEB LOTT
1502	C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
1503	C DU BILAN DE MOMENT ANGULAIRE.
1504	if (bilansmc.eq.1) then
1505	write(,)'Fermeture de aaam_bud.out (FL Vous parle)'
1506	close(27)
1507	close(28)
1508	endif !bilansmc
1509	CMODFIN
1510	c-------------
1511	c--------------SLEBONNOIS
1512	C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
1513	C DES BALLONS
1514	if (ballons.eq.1) then
1515	write(,)'Fermeture des ballons*.out'
1516	close(30)
1517	close(31)
1518	close(32)
1519	close(33)
1520	close(34)
1521	endif !ballons
1522	c-------------
1523	ENDIF
1524
1525
1526	RETURN
1527	END
1528
1529
1530
1531	***********************************************************************
1532	***********************************************************************
1533	***********************************************************************
1534	***********************************************************************
1535	***********************************************************************
1536	***********************************************************************
1537	***********************************************************************
1538	***********************************************************************
1539
1540	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
1541	IMPLICIT none
1542	c
1543	c Tranformer une variable de la grille physique a
1544	c la grille d'ecriture
1545	c
1546	INTEGER nfield,nlon,iim,jjmp1, jjm
1547	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
1548	c
1549	INTEGER i, n, ig
1550	c
1551	jjm = jjmp1 - 1
1552	DO n = 1, nfield
1553	DO i=1,iim
1554	ecrit(i,n) = fi(1,n)
1555	ecrit(i+jjm*iim,n) = fi(nlon,n)
1556	ENDDO
1557	DO ig = 1, nlon - 2
1558	ecrit(iim+ig,n) = fi(1+ig,n)
1559	ENDDO
1560	ENDDO
1561	RETURN
1562	END

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

Download in other formats:

Original Format