Context Navigation

physiq.F @ 97

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

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

Download in other formats:

Original Format