Context Navigation

phyparam.F @ 4200

Last change on this file since 4200 was 4200, checked in by dubos, 5 years ago
simple_physics : cleanup solar angles
File size: 19.9 KB

Line
1	SUBROUTINE phyparam(ngrid,nlayer,nq,
2	s firstcall,lastcall,
3	s rjourvrai,gmtime,ptimestep,
4	s pplev,pplay,pphi,pphis,presnivs,
5	s pu,pv,pt,pq,
6	s pw,
7	s pdu,pdv,pdt,pdq,pdpsrf)
8
9	USE comsaison
10	USE dimphy
11	USE comgeomfi
12	USE phys_const
13	USE planet
14	USE astronomy
15	USE vdif_mod, ONLY : vdif
16	USE solar, ONLY : solang
17	USE radiative, ONLY : mucorr, sw
18	USE radiative_lw, ONLY : lw
19
20	c
21	IMPLICIT NONE
22	c=======================================================================
23	c
24	c subject:
25	c --------
26	c
27	c Organisation of the physical parametrisations of the LMD
28	c 20 parameters GCM for planetary atmospheres.
29	c It includes:
30	c raditive transfer (long and shortwave) for CO2 and dust.
31	c vertical turbulent mixing
32	c convective adjsutment
33	c
34	c author: Frederic Hourdin 15 / 10 /93
35	c -------
36	c
37	c arguments:
38	c ----------
39	c
40	c input:
41	c ------
42	c
43	c ngrid Size of the horizontal grid.
44	c All internal loops are performed on that grid.
45	c nlayer Number of vertical layers.
46	c nq Number of advected fields
47	c firstcall True at the first call
48	c lastcall True at the last call
49	c rjourvrai Number of days counted from the North. Spring
50	c equinoxe.
51	c gmtime hour (s)
52	c ptimestep timestep (s)
53	c pplay(ngrid,nlayer+1) Pressure at the middle of the layers (Pa)
54	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
55	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
56	c pu(ngrid,nlayer) u component of the wind (ms-1)
57	c pv(ngrid,nlayer) v component of the wind (ms-1)
58	c pt(ngrid,nlayer) Temperature (K)
59	c pq(ngrid,nlayer,nq) Advected fields
60	c pudyn(ngrid,nlayer) \
61	c pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
62	c ptdyn(ngrid,nlayer) / corresponding variables
63	c pqdyn(ngrid,nlayer,nq) /
64	c pw(ngrid,?) vertical velocity
65	c
66	c output:
67	c -------
68	c
69	c pdu(ngrid,nlayermx) \
70	c pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding
71	c pdt(ngrid,nlayermx) / variables due to physical processes.
72	c pdq(ngrid,nlayermx) /
73	c pdpsrf(ngrid) /
74	c
75	c=======================================================================
76	c
77	c-----------------------------------------------------------------------
78	c
79	c 0. Declarations :
80	c ------------------
81
82	#include "dimensions.h"
83	#include "callkeys.h"
84	#include "surface.h"
85
86	c Arguments :
87	c -----------
88
89	c inputs:
90	c -------
91	INTEGER ngrid,nlayer,nq
92
93	REAL ptimestep
94	real zdtime
95	REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
96	REAL pphi(ngrid,nlayer)
97	REAL pphis(ngrid)
98	REAL pu(ngrid,nlayer),pv(ngrid,nlayer)
99	REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq)
100	REAL pdu(ngrid,nlayer),pdv(ngrid,nlayer)
101
102	c dynamial tendencies
103	REAL pdtdyn(ngrid,nlayer),pdqdyn(ngrid,nlayer,nq)
104	REAL pdudyn(ngrid,nlayer),pdvdyn(ngrid,nlayer)
105	REAL pw(ngrid,nlayer)
106
107	c Time
108	real rjourvrai
109	REAL gmtime
110
111	c outputs:
112	c --------
113
114	c physical tendencies
115	REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)
116	REAL pdpsrf(ngrid)
117	LOGICAL firstcall,lastcall
118
119	c Local variables :
120	c -----------------
121
122	INTEGER j,l,ig,ierr,aslun,nlevel,igout,it1,it2,isoil,iq
123	INTEGER*4 day_ini
124	c
125	REAL,DIMENSION(ngrid) :: mu0,fract
126	REAL zday
127	REAL zh(ngrid,nlayer),z1,z2
128	REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)
129	REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer)
130	REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid)
131	REAL zflubid(ngrid),zpmer(ngrid)
132	REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
133	REAL zdum1(ngrid,nlayer)
134	REAL zdum2(ngrid,nlayer)
135	REAL zdum3(ngrid,nlayer)
136	REAL ztim1,ztim2,ztim3
137	REAL zls,zinsol
138	REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)
139	REAL zfluxsw(ngrid),zfluxlw(ngrid)
140	character*2 str2
141	REAL factq(nq),tauq(nq)
142	character*3 nomq
143
144	c Local saved variables:
145	c ----------------------
146
147	INTEGER, SAVE :: icount
148	real, SAVE :: zday_last
149	!$OMP THREADPRIVATE( icount,zday_last)
150
151	REAL zps_av
152
153	real,allocatable,SAVE :: tsurf(:),tsoil(:,:),rnatur(:)
154	real,allocatable,SAVE :: capcal(:),fluxgrd(:)
155	real,allocatable,SAVE :: dtrad(:,:),fluxrad(:)
156	real,allocatable,SAVE :: q2(:,:),q2l(:,:)
157	real,allocatable,SAVE :: albedo(:),emissiv(:),z0(:),inertie(:)
158	real,SAVE :: solarcst=1370.
159	real,SAVE :: stephan=5.67e-08
160
161	!$OMP THREADPRIVATE(tsurf,tsoil,rnatur)
162	!$OMP THREADPRIVATE( capcal,fluxgrd,dtrad,fluxrad)
163	!$OMP THREADPRIVATE( q2,q2l)
164	!$OMP THREADPRIVATE( albedo,emissiv,solarcst,z0,inertie)
165	!$OMP THREADPRIVATE( stephan)
166
167
168	EXTERNAL convadj
169	EXTERNAL ismin,ismax
170
171
172	INTEGER longcles
173	PARAMETER ( longcles = 20 )
174	REAL clesphy0( longcles )
175	REAL presnivs(nlayer)
176
177
178
179	print*,'OK DANS PHYPARAM'
180
181	c-----------------------------------------------------------------------
182	c 1. Initialisations :
183	c --------------------
184
185	! call initial0(ngridnlayermxnqmx,pdq)
186	nlevel=nlayer+1
187
188	! print*,'OK ',nlevel
189
190	igout=ngrid/2+1
191	! print*,'OK PHYPARAM ',ngrid,igout
192
193
194	zday=rjourvrai+gmtime
195
196	! print*,'OK PHYPARAM 0A nq ',nq
197	tauq(1)=1800.
198	tauq(2)=10800.
199	tauq(3)=86400.
200	tauq(4)=864000.
201	! print*,'OK PHYPARAM 0 B nq ',nq
202	factq(1:4)=(1.-exp(-ptimestep/tauq(1:4)))/ptimestep
203
204	! print*,'OK PHYPARAM 0 '
205	print*,'nq ',nq
206	print*,'latitude0',ngrid,lati(1:2),lati(ngrid-1:ngrid)
207	print*,'nlayer',nlayer
208	print,'size pdq ',ngridnlayer4,ngridnlayer*nq,
209	s size(pdq),size(lati),size(pq),size(factq)
210	do iq=1,4
211	do l=1,nlayer
212	pdq(1:ngrid,l,iq)=
213	& (1.+sin(lati(1:ngrid))-pq(1:ngrid,l,iq))*factq(iq)
214	enddo
215	enddo
216
217	IF(firstcall) THEN
218
219	print*,'AKk',ngrid,nsoilmx
220	allocate(tsurf(ngrid))
221	print*,'AKa'
222	allocate (tsoil(ngrid,nsoilmx))
223	print*,'AKb'
224	allocate (rnatur(ngrid))
225	print*,'AK2'
226	allocate(capcal(ngrid),fluxgrd(ngrid))
227	print*,'AK3'
228	allocate(dtrad(ngrid,nlayer),fluxrad(ngrid))
229	print*,'AK4'
230	allocate(q2(ngrid,nlayer+1),q2l(ngrid,nlayer+1))
231	print*,'AK5'
232	allocate(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid))
233	print*,'AK6'
234
235
236	do l=1,nlayer
237	pdq(:,l,5)=1.+sin(lati(:))/ptimestep
238	enddo
239	PRINT*,'FIRSTCALL '
240
241	! zday_last=rjourvrai
242	zday_last=zday-ptimestep/unjours
243	c CALL readfi(ngrid,nlayer,nsoilmx,ldrs,
244	c . day_ini,gmtime,albedo,inertie,emissiv,z0,rnatur,
245	c . q2,q2l,tsurf,tsoil)
246	rnatur=1.
247	emissiv(:)=(1.-rnatur(:))emi_mer+rnatur(:)emi_ter
248	inertie(:)=(1.-rnatur(:))I_mer+rnatur(:)I_ter
249	albedo(:)=(1.-rnatur(:))alb_mer+rnatur(:)alb_ter
250	z0(:)=(1.-rnatur(:))Cd_mer+rnatur(:)Cd_ter
251	q2=1.e-10
252	q2l=1.e-10
253	tsurf=300.
254	tsoil=300.
255	print*,tsoil(ngrid/2+1,nsoilmx/2+2)
256	print*,'TS ',tsurf(igout),tsoil(igout,5)
257	CALL iniorbit
258
259	if (.not.callrad) then
260	do ig=1,ngrid
261	fluxrad(ig)=0.
262	enddo
263	endif
264
265	! print*,'OK PHYPARAM 1 '
266	IF(callsoil) THEN
267	CALL soil(ngrid,nsoilmx,firstcall,inertie,
268	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
269	ELSE
270	PRINT*,'WARNING!!! Thermal conduction in the soil
271	s turned off'
272	DO ig=1,ngrid
273	capcal(ig)=1.e5
274	fluxgrd(ig)=0.
275	ENDDO
276	ENDIF
277	c CALL inifrict(ptimestep)
278	icount=0
279
280	PRINT*,'FIRSTCALL B '
281	print*,'INIIO avant iophys_ini '
282	call iophys_ini('phys.nc ',presnivs)
283
284	ENDIF
285	icount=icount+1
286
287	PRINT*,'FIRSTCALL AP '
288	IF (ngrid.NE.ngridmax) THEN
289	PRINT*,'STOP in inifis'
290	PRINT*,'Probleme de dimenesions :'
291	PRINT*,'ngrid = ',ngrid
292	PRINT*,'ngridmax = ',ngridmax
293	STOP
294	ENDIF
295
296	DO l=1,nlayer
297	DO ig=1,ngrid
298	pdv(ig,l)=0.
299	pdu(ig,l)=0.
300	pdt(ig,l)=0.
301	ENDDO
302	ENDDO
303
304	DO ig=1,ngrid
305	pdpsrf(ig)=0.
306	zflubid(ig)=0.
307	zdtsrf(ig)=0.
308	ENDDO
309
310	c-----------------------------------------------------------------------
311	c calcul du geopotentiel aux niveaux intercouches
312	c ponderation des altitudes au niveau des couches en dp/p
313
314	DO l=1,nlayer
315	DO ig=1,ngrid
316	zzlay(ig,l)=pphi(ig,l)/g
317	ENDDO
318	ENDDO
319	DO ig=1,ngrid
320	zzlev(ig,1)=0.
321	ENDDO
322	DO l=2,nlayer
323	DO ig=1,ngrid
324	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
325	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
326	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
327	ENDDO
328	ENDDO
329
330	c-----------------------------------------------------------------------
331	c Transformation de la temperature en temperature potentielle
332	DO l=1,nlayer
333	DO ig=1,ngrid
334	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
335	ENDDO
336	ENDDO
337	DO l=1,nlayer
338	DO ig=1,ngrid
339	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
340	ENDDO
341	ENDDO
342
343	c-----------------------------------------------------------------------
344	c 2. Calcul of the radiative tendencies :
345	c ---------------------------------------
346
347	! print*,'callrad0'
348	IF(callrad) THEN
349	! print*,'callrad'
350
351	c WARNING !!! on calcule le ray a chaque appel
352	c IF( MOD(icount,iradia).EQ.0) THEN
353
354	CALL solarlong(zday,zls)
355	CALL orbite(zls,dist_sol,declin)
356	c declin=0.
357	! print*,'ATTENTIOn : pdeclin = 0',' L_s=',zls
358	print*,'diurnal=',diurnal
359	IF(diurnal) THEN
360	if ( 1.eq.1 ) then
361	ztim1=SIN(declin)
362	ztim2=COS(declin)COS(2.pi*(zday-.5))
363	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
364	c call dump2d(iim,jjm-1,lati(2),'LATI 0 ')
365	c call dump2d(iim,jjm-1,long(2),'LONG 0 ')
366	c call dump2d(iim,jjm-1,sinlon(2),'sinlon0 ')
367	c call dump2d(iim,jjm-1,coslon(2),'coslon0 ')
368	c call dump2d(iim,jjm-1,sinlat(2),'sinlat ')
369	c call dump2d(iim,jjm-1,coslat(2),'coslat ')
370	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
371	s ztim1,ztim2,ztim3,
372	s mu0,fract)
373	else
374	zdtime=ptimestep*float(iradia)
375	CALL zenang(zls,gmtime,zdtime,lati,long,mu0,fract)
376	print*,'ZENANG '
377	endif
378
379	IF(lverbose) THEN
380	PRINT*,'day, declin, sinlon,coslon,sinlat,coslat'
381	PRINT*,zday, declin,
382	s sinlon(igout),coslon(igout),
383	s sinlat(igout),coslat(igout)
384	ENDIF
385	ELSE
386	print*,'declin,ngrid,rad',declin,ngrid,rad
387
388	c call dump2d(iim,jjm-1,lati(2),'LATI ')
389	CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
390	ENDIF
391	c call dump2d(iim,jjm-1,fract(2),'FRACT A ')
392	c call dump2d(iim,jjm-1,mu0(2),'MU0 A ')
393
394
395	c 2.2 Calcul of the radiative tendencies and fluxes:
396	c --------------------------------------------------
397
398	c 2.1.2 levels
399
400	zinsol=solarcst/(dist_sol*dist_sol)
401	print*,iim,jjm,llm,ngrid,nlayer,ngridmax,nlayer
402	print*,'iim,jjm,llm,ngrid,nlayer,ngridmax,nlayer'
403	c call dump2d(iim,jjm-1,albedo(2),'ALBEDO ')
404	c call dump2d(iim,jjm-1,mu0(2),'MU0 ')
405	c call dump2d(iim,jjm-1,fract(2),'FRACT ')
406	c call dump2d(iim,jjm-1,lati(2),'LATI ')
407	zps_av=1.e5
408	if (firstcall) then
409	print*,'WARNING ps_rad impose'
410	endif
411	CALL sw(ngrid,nlayer,diurnal,coefvis,albedo,
412	$ pplev,zps_av,
413	$ mu0,fract,zinsol,
414	$ zfluxsw,zdtsw,
415	$ lverbose)
416	c call dump2d(iim,jjm-1,zfluxsw(2),'SWS 1 ')
417	c stop
418
419	CALL lw(ngrid,nlayer,coefir,emissiv,
420	$ pplev,zps_av,tsurf,pt,
421	$ zfluxlw,zdtlw,
422	$ lverbose)
423
424
425	c 2.4 total flux and tendencies:
426	c ------------------------------
427
428	c 2.4.1 fluxes
429
430	DO ig=1,ngrid
431	fluxrad(ig)=emissiv(ig)*zfluxlw(ig)
432	$ +zfluxsw(ig)*(1.-albedo(ig))
433	zplanck(ig)=tsurf(ig)*tsurf(ig)
434	zplanck(ig)=emissiv(ig)*
435	$ stephanzplanck(ig)zplanck(ig)
436	fluxrad(ig)=fluxrad(ig)-zplanck(ig)
437	ENDDO
438
439	c 2.4.2 temperature tendencies
440
441	DO l=1,nlayer
442	DO ig=1,ngrid
443	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
444	ENDDO
445	ENDDO
446
447	c ENDIF
448
449	c 2.5 Transformation of the radiative tendencies:
450	c -----------------------------------------------
451
452	DO l=1,nlayer
453	DO ig=1,ngrid
454	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
455	ENDDO
456	ENDDO
457
458	IF(lverbose) THEN
459	PRINT*,'Diagnotique for the radaition'
460	PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'
461	PRINT*,albedo(igout),emissiv(igout),mu0(igout),
462	s fract(igout),
463	s fluxrad(igout),zplanck(igout)
464	PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'
465	PRINT*,'unjours',unjours
466	DO l=1,nlayer
467	WRITE(*,'(3f15.5,2e15.2)') pt(igout,l),
468	s pplay(igout,l),pplev(igout,l),
469	s zdtsw(igout,l),zdtlw(igout,l)
470	ENDDO
471	ENDIF
472
473
474	ENDIF
475
476	c-----------------------------------------------------------------------
477	c 3. Vertical diffusion (turbulent mixing):
478	c -----------------------------------------
479	c
480	IF(calldifv) THEN
481
482	DO ig=1,ngrid
483	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
484	ENDDO
485
486	zdum1(:,:)=0.
487	zdum2(:,:)=0.
488
489	do l=1,nlayer
490	do ig=1,ngrid
491	zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l)
492	enddo
493	enddo
494
495	CALL vdif(ngrid,nlayer,zday,
496	$ ptimestep,capcal,z0,
497	$ pplay,pplev,zzlay,zzlev,
498	$ pu,pv,zh,tsurf,emissiv,
499	$ zdum1,zdum2,zdum3,zflubid,
500	$ zdufr,zdvfr,zdhfr,zdtsrfr,q2,q2l,
501	$ lverbose)
502	c igout=ngrid/2+1
503	c PRINT*,'zdufr zdvfr zdhfr'
504	c DO l=1,nlayer
505	c PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l)
506	c ENDDO
507	c CALL difv (ngrid,nlayer,
508	c $ area,lati,capcal,
509	c $ pplev,pphi,
510	c $ pu,pv,zh,tsurf,emissiv,
511	c $ zdum1,zdum2,zdum3,zflubid,
512	c $ zdufr,zdvfr,zdhfr,zdtsrf,
513	c $ .true.)
514	c PRINT*,'zdufr zdvfr zdhfr'
515	c DO l=1,nlayer
516	c PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l)
517	c ENDDO
518	c STOP
519
520	DO l=1,nlayer
521	DO ig=1,ngrid
522	pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
523	pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
524	pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
525	ENDDO
526	ENDDO
527
528	DO ig=1,ngrid
529	zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig)
530	ENDDO
531
532	ELSE
533	DO ig=1,ngrid
534	zdtsrf(ig)=zdtsrf(ig)+
535	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
536	ENDDO
537	ENDIF
538	c
539	c-----------------------------------------------------------------------
540	c 4. Dry convective adjustment:
541	c -----------------------------
542
543	IF(calladj) THEN
544
545	DO l=1,nlayer
546	DO ig=1,ngrid
547	zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l)
548	ENDDO
549	ENDDO
550
551	zdufr(:,:)=0.
552	zdvfr(:,:)=0.
553	zdhfr(:,:)=0.
554
555	CALL convadj(ngrid,nlayer,ptimestep,
556	$ pplay,pplev,zpopsk,
557	$ pu,pv,zh,
558	$ pdu,pdv,zdum1,
559	$ zdufr,zdvfr,zdhfr)
560
561	DO l=1,nlayer
562	DO ig=1,ngrid
563	pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
564	pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
565	pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
566	ENDDO
567	ENDDO
568
569	ENDIF
570
571	c-----------------------------------------------------------------------
572	c On incremente les tendances physiques de la temperature du sol:
573	c ---------------------------------------------------------------
574
575	DO ig=1,ngrid
576	tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)
577	ENDDO
578
579	WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1)
580
581	c-----------------------------------------------------------------------
582	c soil temperatures:
583	c --------------------
584
585	IF (callsoil) THEN
586	CALL soil(ngrid,nsoilmx,.false.,inertie,
587	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
588	IF(lverbose) THEN
589	PRINT*,'Surface Heat capacity,conduction Flux, Ts,
590	s dTs, dt'
591	PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout),
592	s zdtsrf(igout),ptimestep
593	ENDIF
594	ENDIF
595
596	c-----------------------------------------------------------------------
597	c sorties:
598	c --------
599
600	c call dump2d(iim,jjm-1,zfluxsw(2),'SWS 2 ')
601	print*,'zday, zday_last ',zday,zday_last,icount
602	if(abs(zday-zday_last-period_sort)<=ptimestep/unjours/10.) then
603	print*,'zday, zday_last SORTIE ',zday,zday_last
604	zday_last=zday
605	c Ecriture/extension de la coordonnee temps
606
607	do ig=1,ngridmax
608	zpmer(ig)=pplev(ig,1)exp(pphi(ig,1)/(r285.))
609	enddo
610
611	call iophys_ecrit('u',llm,'Vent zonal moy','m/s',pu)
612	call iophys_ecrit('v',llm,'Vent meridien moy','m/s',pv)
613	call iophys_ecrit('temp',llm,'Temperature','K',pt)
614	call iophys_ecrit('geop',llm,'Geopotential','m2/s2',pphi)
615	call iophys_ecrit('plev',llm,'plev','Pa',pplev(:,1:nlayer))
616
617	call iophys_ecrit('du',llm,'du',' ',pdu)
618	call iophys_ecrit('dv',llm,'du',' ',pdv)
619	call iophys_ecrit('dt',llm,'du',' ',pdt)
620	call iophys_ecrit('dtsw',llm,'dtsw',' ',zdtsw)
621	call iophys_ecrit('dtlw',llm,'dtlw',' ',zdtlw)
622
623	do iq=1,nq
624	nomq="tr."
625	write(nomq(2:3),'(i1.1)') iq
626	call iophys_ecrit(nomq,llm,nomq,' ',pq(:,:,iq))
627	enddo
628
629	call iophys_ecrit('ts',1,'Surface temper','K',tsurf)
630	call iophys_ecrit('coslon',1,'coslon',' ',coslon)
631	call iophys_ecrit('sinlon',1,'sinlon',' ',sinlon)
632	call iophys_ecrit('coslat',1,'coslat',' ',coslat)
633	call iophys_ecrit('sinlat',1,'sinlat',' ',sinlat)
634	call iophys_ecrit('mu0',1,'mu0',' ',mu0)
635	call iophys_ecrit('alb',1,'alb',' ',albedo)
636	call iophys_ecrit('fract',1,'fract',' ',fract)
637	call iophys_ecrit('ps',1,'Surface pressure','Pa',pplev)
638	call iophys_ecrit('slp',1,'Sea level pressure','Pa',zpmer)
639	call iophys_ecrit('swsurf',1,'SW surf','Pa',zfluxsw)
640	call iophys_ecrit('lwsurf',1,'LW surf','Pa',zfluxlw)
641
642	endif
643
644	c-----------------------------------------------------------------------
645	IF(lastcall) THEN
646	call iotd_fin
647	PRINT*,'Ecriture du fichier de reinitialiastion de la physique'
648	! if (ierr.ne.0) then
649	! write(6,*)' Pb d''ouverture du fichier restart'
650	! write(6,*)' ierr = ', ierr
651	! call exit(1)
652	! endif
653	write(75) tsurf,tsoil
654	c s (tsurf(1),ig=1,iim+1),
655	c s ( (tsurf(ig),ig=(j-2)iim+2,(j-1)iim+1),
656	c s tsurf((j-2)*iim+2) ,j=2,jjm),
657	c s (tsurf(ngridmax),ig=1,iim+1),
658	c s ( (tsoil(1,l),ig=1,iim+1),
659	c s ( (tsoil(ig,l),ig=(j-2)iim+2,(j-1)iim+1),
660	c s tsoil((j-2)*iim+2,l) ,ig=2,jjm),
661	c s (tsoil(ngridmax,l),ig=1,iim+1)
662	c s ,l=1,nsoilmx)
663	ENDIF
664
665
666	RETURN
667	END

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

Download in other formats:

Original Format