Context Navigation

physiq.F @ 851

Last change on this file since 851 was 849, checked in by emillour, 13 years ago

Work on common dynamics and interfacing with different physics:

Put calls to PVtheta in dynamics between CPP_EARTH flags (because it calls tetalevel, which is supposed to be in the physics; only OK for Earth...).
Adapted makelmdz script so that one can compile main programs in dyn* or phy* (makelmdz_fcm already capable of doing that).
Moved start2archive-VT.F and start2archive-VT.F to phyvenus (as start2archive.F and newstart.F); leave adapting them to Titan for later.
Small correction to phyvenus/testphys1d.F (use module control_mod instead of control.h and call disvert_noterre).
removed "use histcom" in phyvenus/physiq.F and phytitan/physiq.F ; it is not needed since there is already a "use ioipsl" (and it moreover confused makelmdz_fcm...)
Had to add declaration of variable "zlsm1" in phytitan/physiq.F because it is used in "write_hist.h"; but note that it is used while not initialized (but what should it be initialized to?).

File size: 51.5 KB

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

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

Download in other formats:

Original Format