Context Navigation

phyparam.F @ 4196

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

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

Download in other formats:

Original Format