Context Navigation

phyparam.F90 @ 4213

Last change on this file since 4213 was 4213, checked in by dubos, 5 years ago
simple_physics : cleanup phyparam (TBC)
File size: 16.7 KB

Line
1	MODULE phyparam_mod
2	IMPLICIT NONE
3	PRIVATE
4	SAVE
5
6	REAL, PARAMETER :: solarcst=1370., stephan=5.67e-08
7
8	REAL, ALLOCATABLE :: tsurf(:),tsoil(:,:),rnatur(:), &
9	capcal(:),fluxgrd(:), &
10	dtrad(:,:),fluxrad(:), &
11	q2(:,:),q2l(:,:), &
12	albedo(:),emissiv(:),z0(:),inertie(:)
13	!$OMP THREADPRIVATE( tsurf,tsoil,rnatur)
14	!$OMP THREADPRIVATE( capcal,fluxgrd,dtrad,fluxrad)
15	!$OMP THREADPRIVATE( q2,q2l)
16	!$OMP THREADPRIVATE( albedo,emissiv,z0,inertie)
17
18	INTEGER :: icount
19	REAL :: zday_last
20	!$OMP THREADPRIVATE( icount,zday_last)
21
22	PUBLIC :: phyparam
23
24	CONTAINS
25
26	SUBROUTINE phyparam(ngrid,nlayer,nq, &
27	& firstcall,lastcall, &
28	& rjourvrai,gmtime,ptimestep, &
29	& pplev,pplay,pphi,pphis,presnivs, &
30	& pu,pv,pt,pq, &
31	& pw, &
32	& pdu,pdv,pdt,pdq,pdpsrf)
33
34	USE callkeys
35	USE comsaison
36	USE dimphy
37	USE comgeomfi
38	USE phys_const
39	USE planet
40	USE astronomy
41	USE turbulence, ONLY : vdif
42	USE convection, ONLY : convadj
43	USE solar, ONLY : solang, zenang, mucorr
44	USE radiative_sw, ONLY : sw
45	USE radiative_lw, ONLY : lw
46	USE surface
47
48	!=======================================================================
49	! Organisation of the physical parametrisations of the LMD
50	! 20 parameters GCM for planetary atmospheres.
51	! It includes:
52	! raditive transfer (long and shortwave) for CO2 and dust.
53	! vertical turbulent mixing
54	! convective adjsutment
55	!
56	! author: Frederic Hourdin 15 / 10 /93
57	!=======================================================================
58
59	INTEGER, INTENT(IN) :: &
60	ngrid, & ! Size of the horizontal grid.
61	nlayer, & ! Number of vertical layers.
62	nq ! Number of advected fields (tracers)
63	LOGICAL, INTENT(IN) :: &
64	firstcall, & ! True at the first call
65	lastcall ! True at the last call
66	REAL, INTENT(IN) :: &
67	rjourvrai, & ! Number of days counted from the North. Spring equinox
68	gmtime, & ! time of the day in seconds
69	ptimestep, & ! timestep (s)
70	pplev(ngrid,nlayer+1), & ! Pressure at interfaces between layers (pa)
71	pplay(ngrid,nlayer), & ! Pressure at the middle of the layers (Pa)
72	pphi(ngrid,nlayer), & ! Geopotential at the middle of the layers (m2s-2)
73	pphis(ngrid), & ! surface geopotential (unused)
74	pu(ngrid,nlayer), & ! u component of the wind (ms-1)
75	pv(ngrid,nlayer), & ! v component of the wind (ms-1)
76	pw(ngrid,nlayer), & ! vertical velocity (unused)
77	pt(ngrid,nlayer), & ! Temperature (K)
78	pq(ngrid,nlayer,nq) ! Advected fields (unused)
79	REAL, INTENT(OUT) :: & ! output : physical tendencies
80	pdu(ngrid,nlayer), &
81	pdv(ngrid,nlayer), &
82	pdt(ngrid,nlayer), &
83	pdq(ngrid,nlayer,nq), &
84	pdpsrf(ngrid)
85
86	! Local variables :
87	REAL, DIMENSION(ngrid) :: mu0,fract
88	INTEGER :: j,l,ig,ierr,aslun,nlevel,igout,it1,it2,isoil,iq
89	INTEGER*4 day_ini
90	!
91	REAL :: zday, zdtime
92	REAL zh(ngrid,nlayer),z1,z2
93	REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)
94	REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer)
95	REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid)
96	REAL zflubid(ngrid),zpmer(ngrid)
97	REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
98	REAL zdum1(ngrid,nlayer)
99	REAL zdum2(ngrid,nlayer)
100	REAL zdum3(ngrid,nlayer)
101	REAL ztim1,ztim2,ztim3
102	REAL zls,zinsol
103	REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)
104	REAL zfluxsw(ngrid),zfluxlw(ngrid)
105	REAL factq(nq),tauq(nq)
106	character*3 nomq
107
108	! Local saved variables:
109	! ----------------------
110
111	REAL zps_av
112
113	INTEGER longcles
114	PARAMETER ( longcles = 20 )
115	REAL clesphy0( longcles )
116	REAL presnivs(nlayer)
117
118	print*,'OK DANS PHYPARAM'
119	print*,'nq ',nq
120	print*,'latitude0',ngrid,lati(1:2),lati(ngrid-1:ngrid)
121	print*,'nlayer',nlayer
122	print,'size pdq ',ngridnlayer4,ngridnlayer*nq, &
123	& size(pdq),size(lati),size(pq),size(factq)
124
125	nlevel=nlayer+1
126	igout=ngrid/2+1
127	zday=rjourvrai+gmtime
128
129	pdq(:,:,:) = 0. ! we do not use tracers in this physics package
130
131	!-----------------------------------------------------------------------
132	! 1. Initialisations :
133	! --------------------
134
135	IF(firstcall) THEN
136
137	print*,'AKk',ngrid,nsoilmx
138	allocate(tsurf(ngrid))
139	print*,'AKa'
140	allocate (tsoil(ngrid,nsoilmx))
141	print*,'AKb'
142	allocate (rnatur(ngrid))
143	print*,'AK2'
144	allocate(capcal(ngrid),fluxgrd(ngrid))
145	print*,'AK3'
146	allocate(dtrad(ngrid,nlayer),fluxrad(ngrid))
147	print*,'AK4'
148	allocate(q2(ngrid,nlayer+1),q2l(ngrid,nlayer+1))
149	print*,'AK5'
150	allocate(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid))
151	print*,'AK6'
152
153
154	PRINT*,'FIRSTCALL '
155
156	! zday_last=rjourvrai
157	zday_last=zday-ptimestep/unjours
158	rnatur=1.
159	emissiv(:)=(1.-rnatur(:))emi_mer+rnatur(:)emi_ter
160	inertie(:)=(1.-rnatur(:))I_mer+rnatur(:)I_ter
161	albedo(:)=(1.-rnatur(:))alb_mer+rnatur(:)alb_ter
162	z0(:)=(1.-rnatur(:))Cd_mer+rnatur(:)Cd_ter
163	q2=1.e-10
164	q2l=1.e-10
165	tsurf=300.
166	tsoil=300.
167	print*,tsoil(ngrid/2+1,nsoilmx/2+2)
168	print*,'TS ',tsurf(igout),tsoil(igout,5)
169	CALL iniorbit
170
171	if (.not.callrad) then
172	do ig=1,ngrid
173	fluxrad(ig)=0.
174	enddo
175	endif
176
177	! print*,'OK PHYPARAM 1 '
178	IF(callsoil) THEN
179	CALL soil(ngrid,nsoilmx,firstcall,inertie, &
180	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
181	ELSE
182	PRINT*,'WARNING!!! Thermal conduction in the soil turned off'
183	DO ig=1,ngrid
184	capcal(ig)=1.e5
185	fluxgrd(ig)=0.
186	ENDDO
187	ENDIF
188	! CALL inifrict(ptimestep)
189	icount=0
190
191	PRINT*,'FIRSTCALL B '
192	print*,'INIIO avant iophys_ini '
193	call iophys_ini('phys.nc ',presnivs)
194
195	ENDIF
196	icount=icount+1
197
198	PRINT*,'FIRSTCALL AP '
199	IF (ngrid.NE.ngridmax) THEN
200	PRINT*,'STOP in inifis'
201	PRINT*,'Probleme de dimenesions :'
202	PRINT*,'ngrid = ',ngrid
203	PRINT*,'ngridmax = ',ngridmax
204	STOP
205	ENDIF
206
207	DO l=1,nlayer
208	DO ig=1,ngrid
209	pdv(ig,l)=0.
210	pdu(ig,l)=0.
211	pdt(ig,l)=0.
212	ENDDO
213	ENDDO
214
215	DO ig=1,ngrid
216	pdpsrf(ig)=0.
217	zflubid(ig)=0.
218	zdtsrf(ig)=0.
219	ENDDO
220
221	!-----------------------------------------------------------------------
222	! calcul du geopotentiel aux niveaux intercouches
223	! ponderation des altitudes au niveau des couches en dp/p
224
225	DO l=1,nlayer
226	DO ig=1,ngrid
227	zzlay(ig,l)=pphi(ig,l)/g
228	ENDDO
229	ENDDO
230	DO ig=1,ngrid
231	zzlev(ig,1)=0.
232	ENDDO
233	DO l=2,nlayer
234	DO ig=1,ngrid
235	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
236	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
237	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
238	ENDDO
239	ENDDO
240
241	!-----------------------------------------------------------------------
242	! Transformation de la temperature en temperature potentielle
243	DO l=1,nlayer
244	DO ig=1,ngrid
245	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
246	ENDDO
247	ENDDO
248	DO l=1,nlayer
249	DO ig=1,ngrid
250	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
251	ENDDO
252	ENDDO
253
254	!-----------------------------------------------------------------------
255	! 2. Calcul of the radiative tendencies :
256	! ---------------------------------------
257
258	! print*,'callrad0'
259	IF(callrad) THEN
260	! print*,'callrad'
261
262	! WARNING !!! on calcule le ray a chaque appel
263	! IF( MOD(icount,iradia).EQ.0) THEN
264
265	CALL solarlong(zday,zls)
266	CALL orbite(zls,dist_sol,declin)
267	! declin=0.
268	! print*,'ATTENTIOn : pdeclin = 0',' L_s=',zls
269	print*,'diurnal=',diurnal
270	IF(diurnal) THEN
271	if ( 1.eq.1 ) then
272	ztim1=SIN(declin)
273	ztim2=COS(declin)COS(2.pi*(zday-.5))
274	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
275	CALL solang(ngrid,sinlon,coslon,sinlat,coslat, &
276	& ztim1,ztim2,ztim3, &
277	& mu0,fract)
278	else
279	zdtime=ptimestep*float(iradia)
280	CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract)
281	print*,'ZENANG '
282	endif
283
284	IF(lverbose) THEN
285	PRINT*,'day, declin, sinlon,coslon,sinlat,coslat'
286	PRINT*,zday, declin, &
287	& sinlon(igout),coslon(igout), &
288	& sinlat(igout),coslat(igout)
289	ENDIF
290	ELSE
291	print*,'declin,ngrid,rad',declin,ngrid,rad
292	CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
293	ENDIF
294
295	! 2.2 Calcul of the radiative tendencies and fluxes:
296	! --------------------------------------------------
297
298	! 2.1.2 levels
299
300	zinsol=solarcst/(dist_sol*dist_sol)
301	zps_av=1.e5
302	if (firstcall) then
303	print*,'WARNING ps_rad impose'
304	endif
305	CALL sw(ngrid,nlayer,diurnal,coefvis,albedo, &
306	& pplev,zps_av, &
307	& mu0,fract,zinsol, &
308	& zfluxsw,zdtsw, &
309	& lverbose)
310
311	CALL lw(ngrid,nlayer,coefir,emissiv, &
312	& pplev,zps_av,tsurf,pt, &
313	& zfluxlw,zdtlw, &
314	& lverbose)
315
316
317	! 2.4 total flux and tendencies:
318	! ------------------------------
319
320	! 2.4.1 fluxes
321
322	DO ig=1,ngrid
323	fluxrad(ig)=emissiv(ig)*zfluxlw(ig) &
324	& +zfluxsw(ig)*(1.-albedo(ig))
325	zplanck(ig)=tsurf(ig)*tsurf(ig)
326	zplanck(ig)=emissiv(ig)* &
327	& stephanzplanck(ig)zplanck(ig)
328	fluxrad(ig)=fluxrad(ig)-zplanck(ig)
329	ENDDO
330
331	! 2.4.2 temperature tendencies
332
333	DO l=1,nlayer
334	DO ig=1,ngrid
335	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
336	ENDDO
337	ENDDO
338
339	! ENDIF
340
341	! 2.5 Transformation of the radiative tendencies:
342	! -----------------------------------------------
343
344	DO l=1,nlayer
345	DO ig=1,ngrid
346	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
347	ENDDO
348	ENDDO
349
350	IF(lverbose) THEN
351	PRINT*,'Diagnotique for the radaition'
352	PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'
353	PRINT*,albedo(igout),emissiv(igout),mu0(igout), &
354	& fract(igout), &
355	& fluxrad(igout),zplanck(igout)
356	PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'
357	PRINT*,'unjours',unjours
358	DO l=1,nlayer
359	WRITE(*,'(3f15.5,2e15.2)') pt(igout,l), &
360	& pplay(igout,l),pplev(igout,l), &
361	& zdtsw(igout,l),zdtlw(igout,l)
362	ENDDO
363	ENDIF
364
365	ENDIF
366
367	!-----------------------------------------------------------------------
368	! 3. Vertical diffusion (turbulent mixing):
369	! -----------------------------------------
370	!
371	IF(calldifv) THEN
372
373	DO ig=1,ngrid
374	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
375	ENDDO
376
377	zdum1(:,:)=0.
378	zdum2(:,:)=0.
379
380	do l=1,nlayer
381	do ig=1,ngrid
382	zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l)
383	enddo
384	enddo
385
386	CALL vdif(ngrid,nlayer,zday, &
387	& ptimestep,capcal,z0, &
388	& pplay,pplev,zzlay,zzlev, &
389	& pu,pv,zh,tsurf,emissiv, &
390	& zdum1,zdum2,zdum3,zflubid, &
391	& zdufr,zdvfr,zdhfr,zdtsrfr,q2,q2l, &
392	& lverbose)
393
394	DO l=1,nlayer
395	DO ig=1,ngrid
396	pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
397	pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
398	pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
399	ENDDO
400	ENDDO
401
402	DO ig=1,ngrid
403	zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig)
404	ENDDO
405
406	ELSE
407	DO ig=1,ngrid
408	zdtsrf(ig)=zdtsrf(ig)+ &
409	& (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
410	ENDDO
411	ENDIF
412	!
413	!-----------------------------------------------------------------------
414	! 4. Dry convective adjustment:
415	! -----------------------------
416
417	IF(calladj) THEN
418
419	DO l=1,nlayer
420	DO ig=1,ngrid
421	zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l)
422	ENDDO
423	ENDDO
424
425	zdufr(:,:)=0.
426	zdvfr(:,:)=0.
427	zdhfr(:,:)=0.
428
429	CALL convadj(ngrid,nlayer,ptimestep, &
430	& pplay,pplev,zpopsk, &
431	& pu,pv,zh, &
432	& pdu,pdv,zdum1, &
433	& zdufr,zdvfr,zdhfr)
434
435	DO l=1,nlayer
436	DO ig=1,ngrid
437	pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
438	pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
439	pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
440	ENDDO
441	ENDDO
442
443	ENDIF
444
445	!-----------------------------------------------------------------------
446	! On incremente les tendances physiques de la temperature du sol:
447	! ---------------------------------------------------------------
448
449	DO ig=1,ngrid
450	tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)
451	ENDDO
452
453	WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1)
454
455	!-----------------------------------------------------------------------
456	! soil temperatures:
457	! --------------------
458
459	IF (callsoil) THEN
460	CALL soil(ngrid,nsoilmx,.false.,inertie, &
461	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
462	IF(lverbose) THEN
463	PRINT*,'Surface Heat capacity,conduction Flux, Ts, dTs, dt'
464	PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout), &
465	& zdtsrf(igout),ptimestep
466	ENDIF
467	ENDIF
468
469	!-----------------------------------------------------------------------
470	! sorties:
471	! --------
472
473	print*,'zday, zday_last ',zday,zday_last,icount
474	if(abs(zday-zday_last-period_sort)<=ptimestep/unjours/10.) then
475	print*,'zday, zday_last SORTIE ',zday,zday_last
476	zday_last=zday
477	! Ecriture/extension de la coordonnee temps
478
479	do ig=1,ngridmax
480	zpmer(ig)=pplev(ig,1)exp(pphi(ig,1)/(r285.))
481	enddo
482
483	call iophys_ecrit('u',nlayer,'Vent zonal moy','m/s',pu)
484	call iophys_ecrit('v',nlayer,'Vent meridien moy','m/s',pv)
485	call iophys_ecrit('temp',nlayer,'Temperature','K',pt)
486	call iophys_ecrit('geop',nlayer,'Geopotential','m2/s2',pphi)
487	call iophys_ecrit('plev',nlayer,'plev','Pa',pplev(:,1:nlayer))
488
489	call iophys_ecrit('du',nlayer,'du',' ',pdu)
490	call iophys_ecrit('dv',nlayer,'du',' ',pdv)
491	call iophys_ecrit('dt',nlayer,'du',' ',pdt)
492	call iophys_ecrit('dtsw',nlayer,'dtsw',' ',zdtsw)
493	call iophys_ecrit('dtlw',nlayer,'dtlw',' ',zdtlw)
494
495	do iq=1,nq
496	nomq="tr."
497	write(nomq(2:3),'(i1.1)') iq
498	call iophys_ecrit(nomq,nlayer,nomq,' ',pq(:,:,iq))
499	enddo
500
501	call iophys_ecrit('ts',1,'Surface temper','K',tsurf)
502	call iophys_ecrit('coslon',1,'coslon',' ',coslon)
503	call iophys_ecrit('sinlon',1,'sinlon',' ',sinlon)
504	call iophys_ecrit('coslat',1,'coslat',' ',coslat)
505	call iophys_ecrit('sinlat',1,'sinlat',' ',sinlat)
506	call iophys_ecrit('mu0',1,'mu0',' ',mu0)
507	call iophys_ecrit('alb',1,'alb',' ',albedo)
508	call iophys_ecrit('fract',1,'fract',' ',fract)
509	call iophys_ecrit('ps',1,'Surface pressure','Pa',pplev)
510	call iophys_ecrit('slp',1,'Sea level pressure','Pa',zpmer)
511	call iophys_ecrit('swsurf',1,'SW surf','Pa',zfluxsw)
512	call iophys_ecrit('lwsurf',1,'LW surf','Pa',zfluxlw)
513
514	endif
515
516	!-----------------------------------------------------------------------
517	IF(lastcall) THEN
518	call iotd_fin
519	PRINT*,'Ecriture du fichier de reinitialiastion de la physique'
520	write(75) tsurf,tsoil
521	ENDIF
522
523	END SUBROUTINE phyparam
524
525
526	SUBROUTINE alloc_phyparam
527	END SUBROUTINE alloc_phyparam
528	END MODULE phyparam_mod

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

Download in other formats:

Original Format