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

Rev	Line
[4176]	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
[4183]	12	USE phys_const
[4184]	13	USE planet
[4190]	14	USE astronomy
[4188]	15	USE vdif_mod, ONLY : vdif
[4200]	16	USE solar, ONLY : solang
[4198]	17	USE radiative, ONLY : mucorr, sw
[4199]	18	USE radiative_lw, ONLY : lw
	19
[4190]	20	c
[4176]	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
[4188]	168	EXTERNAL convadj
[4176]	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)
[4191]	257	CALL iniorbit
[4176]	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
[4192]	486	zdum1(:,:)=0.
	487	zdum2(:,:)=0.
	488
[4176]	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
[4192]	550
	551	zdufr(:,:)=0.
	552	zdvfr(:,:)=0.
	553	zdhfr(:,:)=0.
	554
[4176]	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