Context Navigation

phyparam.F @ 4181

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

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

Download in other formats:

Original Format