Context Navigation

physiq.F @ 1243

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

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

Download in other formats:

Original Format