Context Navigation

leapfrog_nogcm.F @ 2106

Last change on this file since 2106 was 2106, checked in by emillour, 6 years ago
Add "nogcm" program in dyn3d (only works in serial mode), from ED & FF. EM
File size: 33.4 KB

Line
1	!
2	! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z emillour $
3	!
4	c
5	c
6	SUBROUTINE leapfrog_nogcm(ucov,vcov,teta,ps,masse,phis,q,
7	& time_0)
8
9
10	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
11	#ifdef CPP_IOIPSL
12	use IOIPSL
13	#endif
14	USE infotrac, ONLY: nqtot,ok_iso_verif
15	USE guide_mod, ONLY : guide_main
16	USE write_field, ONLY: writefield
17	USE control_mod, ONLY: planet_type,nday,day_step,iperiod,iphysiq,
18	& less1day,fractday,ndynstep,iconser,
19	& dissip_period,offline,ip_ebil_dyn,
20	& ok_dynzon,periodav,ok_dyn_ave,iecri,
21	& ok_dyn_ins,output_grads_dyn
22	use exner_hyb_m, only: exner_hyb
23	use surfdat_h, only: co2ice
24	use exner_milieu_m, only: exner_milieu
25	use cpdet_mod, only: cpdet,tpot2t,t2tpot
26	use sponge_mod, only: callsponge,mode_sponge,sponge
27	use comuforc_h
28	USE comvert_mod, ONLY: ap,bp,pressure_exner,presnivs,
29	& aps,bps,presnivs,pseudoalt,preff,scaleheight
30	USE comconst_mod, ONLY: daysec,dtvr,dtphys,dtdiss,
31	. cpp,ihf,iflag_top_bound,pi,kappa
32	USE logic_mod, ONLY: iflag_phys,ok_guide,forward,leapf,apphys,
33	. statcl,conser,purmats,tidal,ok_strato
34	USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref,
35	. start_time,dt
36	USE comcstfi_h, ONLY: r,g
37	USE tracer_mod, ONLY: igcm_co2
38
39
40
41
42	IMPLICIT NONE
43
44	c ...... Version du 10/01/98 ..........
45
46	c avec coordonnees verticales hybrides
47	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
48
49	c=======================================================================
50	c
51	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
52	c -------
53	c
54	c Objet:
55	c ------
56	c
57	c GCM LMD nouvelle grille
58	c
59	c=======================================================================
60	c
61	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
62	c et possibilite d'appeler une fonction f(y) a derivee tangente
63	c hyperbolique a la place de la fonction a derivee sinusoidale.
64
65	c ... Possibilite de choisir le shema pour l'advection de
66	c q , en modifiant iadv dans traceur.def (10/02) .
67	c
68	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
69	c Pour Van-Leer iadv=10
70	c
71	c-----------------------------------------------------------------------
72	c Declarations:
73	c -------------
74
75	#include "dimensions.h"
76	#include "paramet.h"
77	#include "comdissnew.h"
78	#include "comgeom.h"
79	!#include "com_io_dyn.h"
80	#include "iniprint.h"
81	#include "academic.h"
82
83	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
84	! #include "clesphys.h"
85
86	REAL,INTENT(IN) :: time_0 ! not used
87
88	c dynamical variables:
89	REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
90	REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
91	REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
92	REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
93	REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
94	REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
95	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
96
97	REAL p (ip1jmp1,llmp1 ) ! interlayer pressure
98	REAL pks(ip1jmp1) ! exner at the surface
99	REAL pk(ip1jmp1,llm) ! exner at mid-layer
100	REAL pkf(ip1jmp1,llm) ! filtered exner at mid-layer
101	REAL phi(ip1jmp1,llm) ! geopotential
102	REAL w(ip1jmp1,llm) ! vertical velocity
103	REAL kpd(ip1jmp1) ! TB15 = exp (-z/H)
104
105	! ADAPTATION GCM POUR CP(T)
106	REAL temp(ip1jmp1,llm) ! temperature
107	REAL tsurpk(ip1jmp1,llm) ! cpp*T/pk
108
109	! real zqmin,zqmax
110
111
112	! ED18 nogcm
113	REAL tau_ps
114	REAL tau_co2
115	REAL tau_teta
116	REAL tetadpmean
117	REAL tetadp(ip1jmp1,llm)
118	REAL dpmean
119	REAL dp(ip1jmp1,llm)
120	REAL tetamean
121	real globaverage2d
122	! LOGICAL firstcall_globaverage2d
123
124
125	c variables dynamiques intermediaire pour le transport
126	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
127
128	c variables dynamiques au pas -1
129	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
130	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
131	REAL massem1(ip1jmp1,llm)
132
133	c tendances dynamiques en */s
134	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
135	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot)
136
137	c tendances de la dissipation en */s
138	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
139	REAL dtetadis(ip1jmp1,llm)
140
141	c tendances de la couche superieure */s
142	c REAL dvtop(ip1jm,llm)
143	REAL dutop(ip1jmp1,llm)
144	c REAL dtetatop(ip1jmp1,llm)
145	c REAL dqtop(ip1jmp1,llm,nqtot),dptop(ip1jmp1)
146
147	c TITAN : tendances due au forces de marees */s
148	REAL dvtidal(ip1jm,llm),dutidal(ip1jmp1,llm)
149
150	c tendances physiques */s
151	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
152	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
153
154	c variables pour le fichier histoire
155	! REAL dtav ! intervalle de temps elementaire
156
157	REAL tppn(iim),tpps(iim),tpn,tps
158	c
159	INTEGER itau,itaufinp1,iav
160	! INTEGER iday ! jour julien
161	REAL time
162
163	REAL SSUM
164	! REAL finvmaold(ip1jmp1,llm)
165
166	cym LOGICAL lafin
167	LOGICAL :: lafin=.false.
168	INTEGER ij,iq,l
169	! INTEGER ik
170
171	! real time_step, t_wrt, t_ops
172
173	REAL rdaym_ini
174	! jD_cur: jour julien courant
175	! jH_cur: heure julienne courante
176	REAL :: jD_cur, jH_cur
177	! INTEGER :: an, mois, jour
178	! REAL :: secondes
179
180	LOGICAL first,callinigrads
181	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
182	save first
183	data first/.true./
184	! real dt_cum
185	! character*10 infile
186	! integer zan, tau0, thoriid
187	! integer nid_ctesGCM
188	! save nid_ctesGCM
189	! real degres
190	! real rlong(iip1), rlatg(jjp1)
191	! real zx_tmp_2d(iip1,jjp1)
192	! integer ndex2d(iip1*jjp1)
193	logical ok_sync
194	parameter (ok_sync = .true.)
195	logical physics
196
197	data callinigrads/.true./
198	character*10 string10
199
200	! REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
201	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
202
203	REAL psmean ! pression moyenne
204	REAL pqco2mean ! moyenne globale ps*qco
205	REAL p0 ! pression de reference
206	REAL p00d ! globalaverage(kpd)
207	REAL qmean_co2,qmean_co2_vert ! mass mean mixing ratio vap co2
208	REAL pqco2(ip1jmp1) ! average co2 mass index : ps*q_co2
209	REAL oldps(ip1jmp1) ! saving old pressure ps to calculate qch4
210
211	c+jld variables test conservation energie
212	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
213	C Tendance de la temp. potentiel d (theta)/ d t due a la
214	C tansformation d'energie cinetique en energie thermique
215	C cree par la dissipation
216	REAL dtetaecdt(ip1jmp1,llm)
217	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
218	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
219	! REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
220	CHARACTER*15 ztit
221	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
222	!IM SAVE ip_ebil_dyn
223	!IM DATA ip_ebil_dyn/0/
224	c-jld
225
226	! integer :: itau_w ! for write_paramLMDZ_dyn.h
227
228	! character*80 dynhist_file, dynhistave_file
229	character(len=*),parameter :: modname="leapfrog"
230	character*80 abort_message
231
232	logical dissip_conservative
233	save dissip_conservative
234	data dissip_conservative/.true./
235
236	INTEGER testita
237	PARAMETER (testita = 9)
238
239	logical , parameter :: flag_verif = .false.
240
241	! for CP(T)
242	real :: dtec
243	real :: ztetaec(ip1jmp1,llm)
244
245	! TEMP : diagnostic mass
246	real :: co2mass(iip1,jjp1)
247	real :: co2ice_ij(iip1,jjp1)
248	integer :: i,j,ig
249	integer, parameter :: ngrid = 2+(jjm-1)*iim
250	! local mass
251	real :: mq(ip1jmp1,llm)
252
253	if (nday>=0) then
254	itaufin = nday*day_step
255	else
256	! to run a given (-nday) number of dynamical steps
257	itaufin = -nday
258	endif
259	if (less1day) then
260	c MODIF VENUS: to run less than one day:
261	itaufin = int(fractday*day_step)
262	endif
263	if (ndynstep.gt.0) then
264	! running a given number of dynamical steps
265	itaufin=ndynstep
266	endif
267	itaufinp1 = itaufin +1
268
269
270
271
272	c INITIALISATIONS
273	dudis(:,:) =0.
274	dvdis(:,:) =0.
275	dtetadis(:,:)=0.
276	dutop(:,:) =0.
277	c dvtop(:,:) =0.
278	c dtetatop(:,:)=0.
279	c dqtop(:,:,:) =0.
280	c dptop(:) =0.
281	dufi(:,:) =0.
282	dvfi(:,:) =0.
283	dtetafi(:,:)=0.
284	dqfi(:,:,:) =0.
285	dpfi(:) =0.
286
287	itau = 0
288	physics=.true.
289	if (iflag_phys==0.or.iflag_phys==2) physics=.false.
290
291	! ED18 nogcm
292	! firstcall_globaverage2d = .true.
293
294	c iday = day_ini+itau/day_step
295	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
296	c IF(time.GT.1.) THEN
297	c time = time-1.
298	c iday = iday+1
299	c ENDIF
300
301
302	c-----------------------------------------------------------------------
303	c On initialise la pression et la fonction d'Exner :
304	c --------------------------------------------------
305
306	dq(:,:,:)=0.
307	CALL pression ( ip1jmp1, ap, bp, ps, p )
308	if (pressure_exner) then
309	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
310	else
311	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
312	endif
313
314	c------------------
315	c TEST PK MONOTONE
316	c------------------
317	write(,) "Test PK"
318	do ij=1,ip1jmp1
319	do l=2,llm
320	! PRINT*,'pk(ij,l) = ',pk(ij,l)
321	! PRINT*,'pk(ij,l-1) = ',pk(ij,l-1)
322	if(pk(ij,l).gt.pk(ij,l-1)) then
323	c write(,) ij,l,pk(ij,l)
324	abort_message = 'PK non strictement decroissante'
325	call abort_gcm(modname,abort_message,1)
326	c write(,) "ATTENTION, Test PK deconnecte..."
327	endif
328	enddo
329	enddo
330	write(,) "Fin Test PK"
331	c stop
332	c------------------
333	c Preparing mixing of pressure and tracers in nogcm
334	c kpd=exp(-z/H) Needed to define a reference pressure :
335	c P0=pmean/globalaverage(kpd)
336	c thus P(i) = P0exp(-z(i)/H) = P0kpd(i)
337	c It is checked that globalaverage2d(Pi)=pmean
338	DO ij=1,ip1jmp1
339	kpd(ij) = exp(-phis(ij)/(r*200.))
340	ENDDO
341	p00d=globaverage2d(kpd) ! mean pres at ref level
342	tau_ps = 1. ! constante de rappel for pressure (s)
343	tau_co2 = 1.E5 !E5 ! constante de rappel for mix ratio qco2 (s)
344	tau_teta = 1.E7 !constante de rappel for potentiel temperature
345
346	! ED18 TEST
347	! PRINT*,'igcm_co2 = ',igcm_co2
348
349	c-----------------------------------------------------------------------
350	c Debut de l'integration temporelle:
351	c ----------------------------------
352
353	1 CONTINUE ! Matsuno Forward step begins here
354
355	c date: (NB: date remains unchanged for Backward step)
356	c -----
357
358	jD_cur = jD_ref + day_ini - day_ref + &
359	& (itau+1)/day_step
360	jH_cur = jH_ref + start_time + &
361	& mod(itau+1,day_step)/float(day_step)
362	jD_cur = jD_cur + int(jH_cur)
363	jH_cur = jH_cur - int(jH_cur)
364
365	c
366
367	! Save fields obtained at previous time step as '...m1'
368	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
369	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
370	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
371	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
372	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
373
374	forward = .TRUE.
375	leapf = .FALSE.
376	dt = dtvr
377
378	2 CONTINUE ! Matsuno backward or leapfrog step begins here
379
380	c-----------------------------------------------------------------------
381
382	c date: (NB: only leapfrog step requires recomputing date)
383	c -----
384
385	IF (leapf) THEN
386	jD_cur = jD_ref + day_ini - day_ref +
387	& (itau+1)/day_step
388	jH_cur = jH_ref + start_time +
389	& mod(itau+1,day_step)/float(day_step)
390	jD_cur = jD_cur + int(jH_cur)
391	jH_cur = jH_cur - int(jH_cur)
392	ENDIF
393
394
395	c gestion des appels de la physique et des dissipations:
396	c ------------------------------------------------------
397	c
398	c ... P.Le Van ( 6/02/95 ) ....
399
400	! ED18: suppression des mentions de la variable apdiss dans le cas
401	! 'nogcm'
402
403	apphys = .FALSE.
404	statcl = .FALSE.
405	conser = .FALSE.
406
407
408	IF( purmats ) THEN
409	! Purely Matsuno time stepping
410	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
411
412
413	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
414	s .and. physics ) apphys = .TRUE.
415	ELSE
416	! Leapfrog/Matsuno time stepping
417	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
418
419
420	IF( MOD(itau+1,iphysiq).EQ.0.AND.physics) apphys=.TRUE.
421	END IF
422
423	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
424	! supress dissipation step
425
426
427
428	c-----------------------------------------------------------------------
429	c calcul des tendances dynamiques:
430	c --------------------------------
431
432	! ED18: suppression de l'onglet pour le cas nogcm
433
434
435	dv(:,:) = 0.
436	du(:,:) = 0.
437	dteta(:,:) = 0.
438	dq(:,:,:) = 0.
439
440
441
442	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
443	c
444	c-----------------------------------------------------------------------
445	c calcul des tendances physiques:
446	c -------------------------------
447	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
448	c
449	IF( purmats ) THEN
450	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
451	ELSE
452	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
453	ENDIF
454	c
455	c
456	IF( apphys ) THEN
457	c
458	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
459	c
460
461	CALL pression ( ip1jmp1, ap, bp, ps, p )
462	if (pressure_exner) then
463	CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf )
464	else
465	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
466	endif
467
468	! Compute geopotential (physics might need it)
469	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
470
471	jD_cur = jD_ref + day_ini - day_ref + &
472	& (itau+1)/day_step
473
474	IF ((planet_type .eq."generic").or.
475	& (planet_type .eq."mars")) THEN
476	! AS: we make jD_cur to be pday
477	jD_cur = int(day_ini + itau/day_step)
478	ENDIF
479
480	! print*,'itau =',itau
481	! print*,'day_step =',day_step
482	! print*,'itau/day_step =',itau/day_step
483
484
485	jH_cur = jH_ref + start_time + &
486	& mod(itau+1,day_step)/float(day_step)
487	IF ((planet_type .eq."generic").or.
488	& (planet_type .eq."mars")) THEN
489	jH_cur = jH_ref + start_time + &
490	& mod(itau,day_step)/float(day_step)
491	ENDIF
492	jD_cur = jD_cur + int(jH_cur)
493	jH_cur = jH_cur - int(jH_cur)
494
495
496
497	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
498	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
499	! write(lunout,*)'current date = ',an, mois, jour, secondes
500
501	c rajout debug
502	c lafin = .true.
503
504
505	c Interface avec les routines de phylmd (phymars ... )
506	c -----------------------------------------------------
507
508	c+jld
509
510	c Diagnostique de conservation de l'Energie : initialisation
511	IF (ip_ebil_dyn.ge.1 ) THEN
512	ztit='bil dyn'
513	! Ehouarn: be careful, diagedyn is Earth-specific!
514	IF (planet_type.eq."earth") THEN
515	CALL diagedyn(ztit,2,1,1,dtphys
516	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
517	ENDIF
518	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
519	c-jld
520	#ifdef CPP_IOIPSL
521	cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
522	cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
523	c IF (first) THEN
524	c first=.false.
525	c#include "ini_paramLMDZ_dyn.h"
526	c ENDIF
527	c
528	c#include "write_paramLMDZ_dyn.h"
529	c
530	#endif
531	! #endif of #ifdef CPP_IOIPSL
532
533	c call WriteField('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
534	! print*,'---------LEAPFROG---------------'
535	! print*,''
536	! print*,'> AVANT CALFIS'
537	! print*,''
538	! print*,'teta(3049,:) = ',teta(3049,:)
539	! print*,''
540	! print*,'dteta(3049,:) = ',dteta(3049,:)
541	! print*,''
542	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
543	! print*,''
544
545
546	CALL calfis( lafin , jD_cur, jH_cur,
547	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
548	$ du,dv,dteta,dq,
549	$ flxw,
550	$ dufi,dvfi,dtetafi,dqfi,dpfi )
551
552	c call WriteField('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
553	c call WriteField('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
554	c call WriteField('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
555	! print*,'> APRES CALFIS (AVANT ADDFI)'
556	! print*,''
557	! print*,'teta(3049,:) = ',teta(3049,:)
558	! print*,''
559	! print*,'dteta(3049,:) = ',dteta(3049,:)
560	! print*,''
561	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
562	! print*,''
563
564
565	c ajout des tendances physiques
566	c ------------------------------
567
568	CALL addfi( dtphys, leapf, forward ,
569	$ ucov, vcov, teta , q ,ps ,
570	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
571
572	! print*,'> APRES ADDFI'
573	! print*,''
574	! print*,'teta(3049,:) = ',teta(3049,:)
575	! print*,''
576	! print*,'dteta(3049,:) = ',dteta(3049,:)
577	! print*,''
578	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
579	! print*,''
580
581	CALL pression (ip1jmp1,ap,bp,ps,p)
582	CALL massdair(p,masse)
583
584	! Mass of tracers in each mess (kg) before Horizontal mixing of pressure
585	c -------------------------------------------------------------------
586
587	DO l=1, llm
588	DO ij=1,ip1jmp1
589	mq(ij,l) = masse(ij,l)*q(ij,l,igcm_co2)
590	ENDDO
591	ENDDO
592
593	c Horizontal mixing of pressure
594	c ------------------------------
595	! Rappel newtonien vers psmean
596	psmean= globaverage2d(ps) ! mean pressure
597	p0=psmean/p00d
598	DO ij=1,ip1jmp1
599	oldps(ij)=ps(ij)
600	ps(ij)=ps(ij) +(p0*kpd(ij)
601	& -ps(ij))*(1-exp(-dtphys/tau_ps))
602	ENDDO
603
604	write(72,*) 'psmean = ',psmean ! mean pressure redistributed
605
606	! security check: test pressure negative
607	DO ij=1,ip1jmp1
608	IF (ps(ij).le.0.) then
609	!write(,) 'Negative pressure :'
610	!write(,) 'ps = ',ps(ij),' ig = ',ij
611	!write(,) 'pmean = ',p0*kpd(ij)
612	!write(,) 'tau_ps = ',tau_ps
613	!STOP
614	ps(ij)=0.0000001*kpd(ij)/p00d
615	ENDIF
616	ENDDO
617	!***********************
618	! Correction on teta due to surface pressure changes
619	DO l = 1,llm
620	DO ij = 1,ip1jmp1
621	teta(ij,l)= teta(ij,l)(ps(ij)/oldps(ij))*kappa
622	ENDDO
623	ENDDO
624	!***********************
625
626
627
628	! ! update pressure and update p and pk
629	! DO ij=1,ip1jmp1
630	! ps(ij) = ps(ij) + dpfi(ij)*dtphys
631	! ENDDO
632	CALL pression (ip1jmp1,ap,bp,ps,p)
633	CALL massdair(p,masse)
634	if (pressure_exner) then
635	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
636	else
637	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
638	endif
639
640	! Update tracers after Horizontal mixing of pressure to ! conserve tracer mass
641	c -------------------------------------------------------------------
642	DO l=1, llm
643	DO ij=1,ip1jmp1
644	q(ij,l,igcm_co2) = mq(ij,l)/ masse(ij,l)
645	ENDDO
646	ENDDO
647
648
649
650	c Horizontal mixing of pressure
651	c ------------------------------
652	! Rappel newtonien vers psmean
653	psmean= globaverage2d(ps) ! mean pressure
654	! ! increment q_co2 with physical tendancy
655	! IF (igcm_co2.ne.0) then
656	! DO l=1, llm
657	! DO ij=1,ip1jmp1
658	! q(ij,l,igcm_co2)=q(ij,l,igcm_co2)+
659	! & dqfi(ij,l,igcm_co2)*dtphys
660	! ENDDO
661	! ENDDO
662	! ENDIF
663
664	c Mixing CO2 vertically
665	c --------------------------
666	if (igcm_co2.ne.0) then
667	DO ij=1,ip1jmp1
668	qmean_co2_vert=0.
669	DO l=1, llm
670	qmean_co2_vert= qmean_co2_vert
671	& + q(ij,l,igcm_co2)*( p(ij,l) - p(ij,l+1))
672	END DO
673	qmean_co2_vert= qmean_co2_vert/ps(ij)
674	DO l=1, llm
675	q(ij,l,igcm_co2)= qmean_co2_vert
676	END DO
677	END DO
678	end if
679
680
681
682	c Horizontal mixing of pressure
683	c ------------------------------
684	! Rappel newtonien vers psmean
685	psmean= globaverage2d(ps) ! mean pressure
686
687	c Horizontal mixing tracers, and temperature (replace dynamics in nogcm)
688	c ---------------------------------------------------------------------------------
689
690	! Simulate redistribution by dynamics for qco2
691	if (igcm_co2.ne.0) then
692
693	DO ij=1,ip1jmp1
694	pqco2(ij)= ps(ij) * q(ij,1,igcm_co2)
695	ENDDO
696	pqco2mean=globaverage2d(pqco2)
697
698	! Rappel newtonien vers qco2_mean
699	qmean_co2= pqco2mean / psmean
700
701	DO ij=1,ip1jmp1
702	q(ij,1,igcm_co2)=q(ij,1,igcm_co2)+
703	& (qmean_co2-q(ij,1,igcm_co2))*
704	& (1.-exp(-dtphys/tau_co2))
705	ENDDO
706
707	DO l=2, llm
708	DO ij=1,ip1jmp1
709	q(ij,l,igcm_co2)=q(ij,1,igcm_co2)
710	END DO
711	END DO
712
713	! TEMPORAIRE (ED)
714	! PRINT*,'psmean = ',psmean
715	! PRINT*,'qmean_co2 = ',qmean_co2
716	! PRINT*,'pqco2mean = ',pqco2mean
717	! PRINT*,'q(50,1,igcm_co2) = ',q(50,1,igcm_co2)
718	! PRINT*,'q(50,2,igcm_co2) = ',q(50,2,igcm_co2)
719	! PRINT*,'q(50,3,igcm_co2) = ',q(50,3,igcm_co2)
720
721	endif ! igcm_co2.ne.0
722
723
724	! ********************************************************
725
726	c Horizontal mixing of pressure
727	c ------------------------------
728	! Rappel newtonien vers psmean
729	psmean= globaverage2d(ps) ! mean pressure
730
731
732	c Mixing Temperature horizontally
733	c -------------------------------
734	! initialize variables that will be averaged
735	DO l=1,llm
736	DO ij=1,ip1jmp1
737	dp(ij,l) = p(ij,l) - p(ij,l+1)
738	tetadp(ij,l) = teta(ij,l)*dp(ij,l)
739	ENDDO
740	ENDDO
741
742	DO l=1,llm
743	tetadpmean = globaverage2d(tetadp(:,l))
744	dpmean = globaverage2d(dp(:,l))
745	tetamean = tetadpmean / dpmean
746	DO ij=1,ip1jmp1
747	teta(ij,l) = teta(ij,l) + (tetamean - teta(ij,l)) *
748	& (1 - exp(-dtphys/tau_teta))
749	ENDDO
750	ENDDO
751
752
753	ENDIF ! of IF( apphys )
754
755
756	c ********************************************************************
757	c ********************************************************************
758	c .... fin de l'integration dynamique et physique pour le pas itau ..
759	c ********************************************************************
760	c ********************************************************************
761
762
763	IF ( .NOT.purmats ) THEN
764	c ........................................................
765	c .............. schema matsuno + leapfrog ..............
766	c ........................................................
767
768	IF(forward. OR. leapf) THEN
769	itau= itau + 1
770	c iday= day_ini+itau/day_step
771	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
772	c IF(time.GT.1.) THEN
773	c time = time-1.
774	c iday = iday+1
775	c ENDIF
776	ENDIF
777
778
779	IF( itau. EQ. itaufinp1 ) then
780	if (flag_verif) then
781	write(79,*) 'ucov',ucov
782	write(80,*) 'vcov',vcov
783	write(81,*) 'teta',teta
784	write(82,*) 'ps',ps
785	write(83,*) 'q',q
786	WRITE(85,*) 'q1 = ',q(:,:,1)
787	WRITE(86,*) 'q3 = ',q(:,:,3)
788	endif
789
790	abort_message = 'Simulation finished'
791
792	call abort_gcm(modname,abort_message,0)
793	ENDIF
794	c-----------------------------------------------------------------------
795	c ecriture du fichier histoire moyenne:
796	c -------------------------------------
797
798	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
799	IF(itau.EQ.itaufin) THEN
800	iav=1
801	ELSE
802	iav=0
803	ENDIF
804
805	! ! Ehouarn: re-compute geopotential for outputs
806	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
807
808	IF (ok_dynzon) THEN
809	#ifdef CPP_IOIPSL
810	c les traceurs ne sont pas sortis, trop lourd.
811	c Peut changer eventuellement si besoin.
812	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
813	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
814	& du,dudis,dutop,dufi)
815	#endif
816	END IF
817	IF (ok_dyn_ave) THEN
818	#ifdef CPP_IOIPSL
819	CALL writedynav(itau,vcov,
820	& ucov,teta,pk,phi,q,masse,ps,phis)
821	#endif
822	ENDIF
823
824	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
825
826	if (ok_iso_verif) then
827	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
828	endif !if (ok_iso_verif) then
829
830	c-----------------------------------------------------------------------
831	c ecriture de la bande histoire:
832	c ------------------------------
833
834	IF( MOD(itau,iecri).EQ.0) THEN
835	! Ehouarn: output only during LF or Backward Matsuno
836	if (leapf.or.(.not.leapf.and.(.not.forward))) then
837	! ADAPTATION GCM POUR CP(T)
838	call tpot2t(ijp1llm,teta,temp,pk)
839	tsurpk = cpp*temp/pk
840	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
841	unat=0.
842	do l=1,llm
843	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
844	vnat(:,l)=vcov(:,l)/cv(:)
845	enddo
846	#ifdef CPP_IOIPSL
847	if (ok_dyn_ins) then
848	! write(lunout,*) "leapfrog: call writehist, itau=",itau
849	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
850	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
851	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
852	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
853	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
854	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
855	endif ! of if (ok_dyn_ins)
856	#endif
857	! For some Grads outputs of fields
858	if (output_grads_dyn) then
859	#include "write_grads_dyn.h"
860	endif
861	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
862	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
863
864	IF(itau.EQ.itaufin) THEN
865
866	if (planet_type=="mars") then
867	CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
868	& vcov,ucov,teta,q,masse,ps)
869	else
870	CALL dynredem1("restart.nc",start_time,
871	& vcov,ucov,teta,q,masse,ps)
872	endif
873	CLOSE(99)
874	!!! Ehouarn: Why not stop here and now?
875	ENDIF ! of IF (itau.EQ.itaufin)
876
877	c-----------------------------------------------------------------------
878	c gestion de l'integration temporelle:
879	c ------------------------------------
880
881	IF( MOD(itau,iperiod).EQ.0 ) THEN
882	GO TO 1
883	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
884
885	IF( forward ) THEN
886	c fin du pas forward et debut du pas backward
887
888	forward = .FALSE.
889	leapf = .FALSE.
890	GO TO 2
891
892	ELSE
893	c fin du pas backward et debut du premier pas leapfrog
894
895	leapf = .TRUE.
896	dt = 2.*dtvr
897	GO TO 2
898	END IF ! of IF (forward)
899	ELSE
900
901	c ...... pas leapfrog .....
902
903	leapf = .TRUE.
904	dt = 2.*dtvr
905	GO TO 2
906	END IF ! of IF (MOD(itau,iperiod).EQ.0)
907	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
908
909	ELSE ! of IF (.not.purmats)
910
911	if (ok_iso_verif) then
912	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
913	endif !if (ok_iso_verif) then
914
915	c ........................................................
916	c .............. schema matsuno ...............
917	c ........................................................
918	IF( forward ) THEN
919
920	itau = itau + 1
921	c iday = day_ini+itau/day_step
922	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
923	c
924	c IF(time.GT.1.) THEN
925	c time = time-1.
926	c iday = iday+1
927	c ENDIF
928
929	forward = .FALSE.
930	IF( itau. EQ. itaufinp1 ) then
931	abort_message = 'Simulation finished'
932	call abort_gcm(modname,abort_message,0)
933	ENDIF
934	GO TO 2
935
936	ELSE ! of IF(forward) i.e. backward step
937
938	if (ok_iso_verif) then
939	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
940	endif !if (ok_iso_verif) then
941
942	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
943	IF(itau.EQ.itaufin) THEN
944	iav=1
945	ELSE
946	iav=0
947	ENDIF
948
949	! ! Ehouarn: re-compute geopotential for outputs
950	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
951
952	IF (ok_dynzon) THEN
953	#ifdef CPP_IOIPSL
954	c les traceurs ne sont pas sortis, trop lourd.
955	c Peut changer eventuellement si besoin.
956	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
957	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
958	& du,dudis,dutop,dufi)
959	#endif
960	ENDIF
961	IF (ok_dyn_ave) THEN
962	#ifdef CPP_IOIPSL
963	CALL writedynav(itau,vcov,
964	& ucov,teta,pk,phi,q,masse,ps,phis)
965	#endif
966	ENDIF
967
968	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
969
970	IF(MOD(itau,iecri ).EQ.0) THEN
971	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
972	! ADAPTATION GCM POUR CP(T)
973	call tpot2t(ijp1llm,teta,temp,pk)
974	tsurpk = cpp*temp/pk
975	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
976	unat=0.
977	do l=1,llm
978	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
979	vnat(:,l)=vcov(:,l)/cv(:)
980	enddo
981	#ifdef CPP_IOIPSL
982	if (ok_dyn_ins) then
983	! write(lunout,*) "leapfrog: call writehist (b)",
984	! & itau,iecri
985	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
986	endif ! of if (ok_dyn_ins)
987	#endif
988	! For some Grads outputs
989	if (output_grads_dyn) then
990	#include "write_grads_dyn.h"
991	endif
992
993	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
994
995	IF(itau.EQ.itaufin) THEN
996	if (planet_type=="mars") then
997	CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
998	& vcov,ucov,teta,q,masse,ps)
999	else
1000	CALL dynredem1("restart.nc",start_time,
1001	& vcov,ucov,teta,q,masse,ps)
1002	endif
1003	ENDIF ! of IF(itau.EQ.itaufin)
1004
1005	forward = .TRUE.
1006	GO TO 1
1007
1008	ENDIF ! of IF (forward)
1009
1010	END IF ! of IF(.not.purmats)
1011
1012	STOP
1013	END
1014
1015	! ************************************************************
1016	! globalaverage VERSION PLuto
1017	! ************************************************************
1018	FUNCTION globaverage2d(var)
1019	! ************************************************************
1020	IMPLICIT NONE
1021	#include "dimensions.h"
1022	#include "paramet.h"
1023	#include "comgeom2.h"
1024	!ip1jmp1 called in paramet.h = iip1 x jjp1
1025	REAL var(iip1,jjp1), globaverage2d
1026	integer i,j
1027	REAL airetot
1028	save airetot
1029	logical firstcall
1030	data firstcall/.true./
1031	save firstcall
1032
1033	if (firstcall) then
1034	firstcall=.false.
1035	airetot =0.
1036	DO j=2,jjp1-1 ! lat
1037	DO i = 1, iim ! lon
1038	airetot = airetot + aire(i,j)
1039	END DO
1040	END DO
1041	DO i=1,iim
1042	airetot=airetot + aire(i,1)+aire(i,jjp1)
1043	ENDDO
1044	end if
1045
1046	globaverage2d = 0.
1047	DO j=2,jjp1-1
1048	DO i = 1, iim
1049	globaverage2d = globaverage2d + var(i,j)*aire(i,j)
1050	END DO
1051	END DO
1052
1053	DO i=1,iim
1054	globaverage2d = globaverage2d + var(i,1)*aire(i,1)
1055	globaverage2d = globaverage2d + var(i,jjp1)*aire(i,jjp1)
1056	ENDDO
1057
1058	globaverage2d = globaverage2d/airetot
1059	return
1060	end
1061

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

Download in other formats:

Original Format