Context Navigation

source: codes/icosagcm/trunk/src/phyparam.F @ 163

Last change on this file since 163 was 163, checked in by dubos, 11 years ago
Disable zenang to enable compilation with -debug
File size: 17.8 KB

Line
1	MODULE PHY
2	USE dimphys_mod
3	LOGICAL:: firstcall,lastcall
4	contains
5
6	SUBROUTINE phyparam_lmd(it,ngrid,nlayer,nq,
7	s ptimestep,lati,long,rjourvrai,gmtime,
8	s pplev,pplay,pphi,pphis,
9	s pu,pv,pt,pq,
10	s pdu,pdv,pdt,pdq,pdpsrf)
11
12	USE ICOSA
13	USE dimphys_mod
14	USE RADIATION
15	USE SURFACE_PROCESS
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,nlayer) \
66	c pdv(ngrid,nlayer) \ Temporal derivative of the corresponding
67	c pdt(ngrid,nlayer) / variables due to physical processes.
68	c pdq(ngrid,nlayer) /
69	c pdpsrf(ngrid) /
70	c
71	c=======================================================================
72	c
73	!c-----------------------------------------------------------------------
74	!c
75	!c 0. Declarations :
76	!c ------------------
77	!c Arguments :
78	!c -----------
79
80	!c inputs:
81	!c -------
82	INTEGER ngrid,nlayer,nq,it,ij,i
83	REAL ptimestep
84	real zdtime
85	REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
86	REAL pphi(ngrid,nlayer)
87	REAL pphis(ngrid)
88	REAL pu(ngrid,nlayer),pv(ngrid,nlayer)
89	REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq)
90	REAL pdu(ngrid,nlayer),pdv(ngrid,nlayer)
91
92	!c dynamial tendencies
93	REAL pdtdyn(ngrid,nlayer),pdqdyn(ngrid,nlayer,nq)
94	REAL pdudyn(ngrid,nlayer),pdvdyn(ngrid,nlayer)
95	REAL pw(ngrid,nlayer)
96
97	!c Time
98	real rjourvrai
99	REAL gmtime
100
101	!c outputs:
102	!c --------
103
104	!c physical tendencies
105	REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)
106	REAL pdpsrf(ngrid)
107
108
109	!c Local variables :
110	!c -----------------
111
112	INTEGER j,l,ig,ierr,aslun,nlevel,igout,it1,it2,unit,isoil
113	REAL zps_av
114	REAL zday
115	REAL zh(ngrid,nlayer),z1,z2
116	REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)
117	REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer)
118	REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid)
119	REAL zflubid(ngrid),zpmer(ngrid)
120	REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
121	REAL zdum1(ngrid,nlayer)
122	REAL zdum2(ngrid,nlayer)
123	REAL zdum3(ngrid,nlayer)
124	REAL ztim1,ztim2,ztim3
125	REAL zls,zinsol
126	REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)
127	REAL zfluxsw(ngrid),zfluxlw(ngrid)
128	character*2 str2
129
130	!c Local saved variables:
131	!c ----------------------
132	REAL(rstd)::long(ngrid),lati(ngrid)
133	REAL(rstd)::mu0(ngrid),fract(ngrid),coslat(ngrid)
134	REAL(rstd)::sinlon(ngrid),coslon(ngrid),sinlat(ngrid)
135	REAL(rstd)::dist_sol,declin
136	REAL::totarea !sarvesh
137	!!!!!!!!sarvesh !!!!!!! CHECK SAVE ATTRIBUTE
138	INTEGER:: icount
139	real:: zday_last
140	REAL:: solarcst
141
142	SAVE icount,zday_last
143	SAVE solarcst
144	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
145	REAL stephan
146	SAVE stephan
147	DATA stephan/5.67e-08/
148	DATA solarcst/1370./
149	REAL presnivs(nlayer)
150	INTEGER:: nn1,nn2
151	c-----------------------------------------------------------------------
152	c 1. Initialisations :
153	c --------------------
154	c call initial0(ngridnlayernqmx,pdq)
155
156	! nn1=(jj_begin -1)*iim+ii_begin
157	! nn2=(jj_end -1)*iim+ii_end
158
159	nlevel=nlayer+1
160	igout=ngrid/2+1
161
162	DO j=jj_begin-offset,jj_end+offset
163	DO i=ii_begin-offset,ii_end+offset
164	ig=(j-1)*iim+i
165	sinlat(ig) = sin(lati(ig))
166	coslat(ig) = cos(lati(ig))
167	sinlon(ig) = sin(long(ig))
168	coslon(ig) = cos(long(ig))
169	END DO
170	ENDDO
171	zday=rjourvrai+gmtime
172	IF ( it == 0 ) then
173	firstcall=.TRUE.
174	ELSE
175	firstcall=.FALSE.
176	ENDIF
177	IF ( it == ndays*day_step ) Then
178	lastcall = .True.
179	END IF
180
181	IF(firstcall) THEN
182	PRINT*,'FIRSTCALL ',ngridmx,nlayermx,nsoilmx
183	zday_last=rjourvrai
184	inertie=2000
185	albedo=0.2
186	emissiv=1.
187	z0=0.1
188	rnatur=1.
189	q2=1.e-10
190	q2l=1.e-10
191	tsurf(:)=300.
192	tsoil(:,:)=300.
193	! print*,tsoil(ngrid/2+1,nsoilmx/2+2)
194	! print*,'TS ',tsurf(igout),tsoil(igout,5)
195
196	IF (.not.callrad) then
197	DO j=jj_begin-offset,jj_end+offset
198	DO i=ii_begin-offset,ii_end+offset
199	ig=(j-1)*iim+i
200	fluxrad(ig)=0.
201	enddo
202	enddo
203	ENDIF
204
205	IF(callsoil) THEN
206	CALL soil(ngrid,nsoilmx,firstcall,inertie,
207	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
208	ELSE
209	PRINT*,'WARNING!!! Thermal conduction in the soil
210	s turned off'
211	DO j=jj_begin-offset,jj_end+offset
212	DO i=ii_begin-offset,ii_end+offset
213	ig=(j-1)*iim+i
214	capcal(ig)=1.e5
215	fluxgrd(ig)=0.
216	ENDDO
217	ENDDO
218	ENDIF
219	icount=0
220	ENDIF
221
222	IF (ngrid.NE.ngrid) THEN
223	PRINT*,'STOP in inifis'
224	PRINT*,'Probleme de dimenesions :'
225	PRINT*,'ngrid = ',ngrid
226	PRINT*,'ngrid = ',ngrid
227	STOP
228	ENDIF
229
230	DO l=1,nlayer
231	DO j=jj_begin-offset,jj_end+offset
232	DO i=ii_begin-offset,ii_end+offset
233	ig=(j-1)*iim+i
234	pdv(ig,l)=0.
235	pdu(ig,l)=0.
236	pdt(ig,l)=0.
237	ENDDO
238	ENDDO
239	ENDDO
240
241	DO j=jj_begin-offset,jj_end+offset
242	DO i=ii_begin-offset,ii_end+offset
243	ig=(j-1)*iim+i
244	pdpsrf(ig)=0.
245	zflubid(ig)=0.
246	zdtsrf(ig)=0.
247	ENDDO
248	ENDDO
249
250	zps_av=0.0
251	DO j=jj_begin-offset,jj_end+offset
252	DO i=ii_begin-offset,ii_end+offset
253	ig=(j-1)*iim+i
254	zps_av=zps_av+pplev(ig,1)*Ai(ig)
255	totarea=totarea+Ai(ig)
256	END DO
257	END DO
258	zps_av=zps_av/totarea
259
260
261	!print*,"maxplev",maxval(pplev(:,1)),minval(pplev(:,1))
262	c-----------------------------------------------------------------------
263	c calcul du geopotentiel aux niveaux intercouches
264	c ponderation des altitudes au niveau des couches en dp/p
265
266	DO l=1,nlayer
267	DO j=jj_begin-offset,jj_end+offset
268	DO i=ii_begin-offset,ii_end+offset
269	ig=(j-1)*iim+i
270	zzlay(ig,l)=pphi(ig,l)/g
271	ENDDO
272	ENDDO
273	ENDDO
274
275	!print*,"zzlay",maxval(zzlay(:,1)),minval(zzlay(:,1))
276
277
278	DO j=jj_begin-offset,jj_end+offset
279	DO i=ii_begin-offset,ii_end+offset
280	ig=(j-1)*iim+i
281	zzlev(ig,1)=0.
282	ENDDO
283	ENDDO
284
285	DO l=2,nlayer
286	DO j=jj_begin-offset,jj_end+offset
287	DO i=ii_begin-offset,ii_end+offset
288	ig=(j-1)*iim+i
289	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
290	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
291	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
292	ENDDO
293	ENDDO
294	ENDDO
295	!print*,"zzlev",maxval(zzlev(:,1)),minval(zzlev(:,1))
296	c-----------------------------------------------------------------------
297	c Transformation de la temperature en temperature potentielle
298	DO l=1,nlayer
299	DO j=jj_begin-offset,jj_end+offset
300	DO i=ii_begin-offset,ii_end+offset
301	ig=(j-1)*iim +i
302	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**kappa
303	ENDDO
304	ENDDO
305	ENDDO
306
307	DO l=1,nlayer
308	DO j=jj_begin-offset,jj_end+offset
309	DO i=ii_begin-offset,ii_end+offset
310	ig=(j-1)*iim+i
311	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
312	ENDDO
313	ENDDO
314	ENDDO
315	!print*,"ph pot",maxval(zh(:,1)),minval(zh(:,1))
316	! go to 101
317	c-----------------------------------------------------------------------
318	c 2. Calcul of the radiative tendencies :
319	c ---------------------------------------
320	! print*,'callrad0'
321	IF(callrad) THEN
322	! print*,'callrad'
323	! WARNING !!! on calcule le ray a chaque appel
324	! IF( MOD(icount,iradia).EQ.0) THEN
325
326	CALL solarlong(zday,zls)
327	CALL orbite(zls,dist_sol,declin)
328	!
329	! print*,'ATTENTIOn : pdeclin = 0',' L_s=',zls
330	! print*,'diurnal=',diurnal
331	IF(diurnal) THEN
332	if ( 1.eq.1 ) then
333	ztim1=SIN(declin)
334	ztim2=COS(declin)COS(2.pi*(zday-.5))
335	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
336
337	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
338	s ztim1,ztim2,ztim3,
339	s mu0,fract)
340	else
341	zdtime=ptimestep*float(iradia)
342	! CALL zenang(zls,gmtime,zdtime,lati,long,mu0,fract) ! FIXME
343	!print*,'ZENANG '
344	endif
345
346	IF(lverbose) THEN
347	PRINT*,'day, declin, sinlon,coslon,sinlat,coslat'
348	PRINT*,zday, declin,
349	s sinlon(igout),coslon(igout),
350	s sinlat(igout),coslat(igout)
351	ENDIF
352	ELSE
353	!print*,'declin,ngrid,radius',declin,ngrid,radius
354	CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,radius)
355	! open(100,file="mu0.txt")
356	! write(100,*)(mu0(ij),ij=1,ngrid)
357	ENDIF
358	! print*,"iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii"
359	!c 2.2 Calcul of the radiative tendencies and fluxes:
360	!c --------------------------------------------------
361	!c 2.1.2 levels
362
363	zinsol=solarcst/(dist_sol*dist_sol)
364	! print*,'iim,jjm,llm,ngrid,nlayer,ngrid,nlayer'
365	! print*,iim,jjm,llm,ngrid,nlayer,ngrid,nlayer
366	! print*,"zinsol sol_dist",zinsol,dist_sol
367	! STOP
368
369	CALL sw(ngrid,nlayer,diurnal,coefvis,albedo,
370	$ pplev,zps_av,
371	$ mu0,fract,zinsol,
372	$ zfluxsw,zdtsw,
373	$ lverbose)
374
375	!print*,"sw",maxval(zfluxsw),minval(zfluxsw),
376	! $ maxval(zdtsw),minval(zdtsw), " it",it
377
378	! print*,"lllllllllllllllllllllllllllllllllllllllll"
379	! print*,"pplev",maxval(pplev),minval(pplev)
380	! print*,"zps,tsurf",zps_av,maxval(tsurf),minval(tsurf)
381	! print*,"pt",maxval(pt),minval(pt)
382
383
384	CALL lw(ngrid,nlayer,coefir,emissiv,
385	$ pplev,zps_av,tsurf,pt,
386	$ zfluxlw,zdtlw,
387	$ lverbose)
388
389	!print*,"lw",maxval(zfluxlw),minval(zfluxlw),
390	! $ maxval(zdtlw),minval(zdtlw)
391
392	! print*,"lw",maxval(
393
394	! print*,"after lw"
395	c 2.4 total flux and tendencies:
396	c ------------------------------
397
398	c 2.4.1 fluxes
399
400	DO j=jj_begin-offset,jj_end+offset
401	DO i=ii_begin-offset,ii_end+offset
402	ig=(j-1)*iim+i
403	fluxrad(ig)=emissiv(ig)*zfluxlw(ig)
404	$ +zfluxsw(ig)*(1.-albedo(ig))
405
406	zplanck(ig)=tsurf(ig)*tsurf(ig)
407
408	zplanck(ig)=emissiv(ig)*
409	$ stephanzplanck(ig)zplanck(ig)
410
411	fluxrad(ig)=fluxrad(ig)-zplanck(ig)
412	ENDDO
413	ENDDO
414	c 2.4.2 temperature tendencies
415
416	DO l=1,nlayer
417	DO j=jj_begin-offset,jj_end+offset
418	DO i=ii_begin-offset,ii_end+offset
419	ig=(j-1)*iim+i
420	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
421	ENDDO
422	ENDDO
423	ENDDO
424
425
426	c 2.5 Transformation of the radiative tendencies:
427	c -----------------------------------------------
428
429	DO l=1,nlayer
430	DO j=jj_begin-offset,jj_end+offset
431	DO i=ii_begin-offset,ii_end+offset
432	ig=(j-1)*iim+i
433	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
434	ENDDO
435	ENDDO
436	ENDDO
437
438	IF(lverbose) THEN
439	PRINT*,'Diagnotique for the radaition'
440	PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'
441	PRINT*,albedo(igout),emissiv(igout),mu0(igout),
442	s fract(igout),
443	s fluxrad(igout),zplanck(igout)
444	PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'
445	! PRINT*,'unjours',unjours
446	DO l=1,nlayer
447	WRITE(*,'(3f15.5,2e15.2)') pt(igout,l),
448	s pplay(igout,l),pplev(igout,l),
449	s zdtsw(igout,l),zdtlw(igout,l)
450	ENDDO
451	ENDIF
452	ENDIF !( CALL RADIATION )
453	! print*,"eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee"
454
455	c-----------------------------------------------------------------------
456	c 3. Vertical diffusion (turbulent mixing):
457	c -----------------------------------------
458	c
459	IF(calldifv) THEN
460
461	DO ig=1,ngrid
462	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
463	ENDDO
464
465	CALL zerophys(ngrid*nlayer,zdum1)
466	CALL zerophys(ngrid*nlayer,zdum2)
467	c CALL zerophys(ngrid*nlayer,zdum3)
468	do l=1,nlayer
469	do ig=1,ngrid
470	zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l)
471	enddo
472	enddo
473
474	CALL vdif(ngrid,nlayer,zday,
475	$ ptimestep,capcal,z0,
476	$ pplay,pplev,zzlay,zzlev,
477	$ pu,pv,zh,tsurf,emissiv,
478	$ zdum1,zdum2,zdum3,zflubid,
479	$ zdufr,zdvfr,zdhfr,zdtsrfr,q2,q2l,
480	$ lverbose)
481	c igout=ngrid/2+1
482	c PRINT*,'zdufr zdvfr zdhfr'
483	c DO l=1,nlayer
484	c PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l)
485	c ENDDO
486	c CALL difv (ngrid,nlayer,
487	c $ area,lati,capcal,
488	c $ pplev,pphi,
489	c $ pu,pv,zh,tsurf,emissiv,
490	c $ zdum1,zdum2,zdum3,zflubid,
491	c $ zdufr,zdvfr,zdhfr,zdtsrf,
492	c $ .true.)
493	c PRINT*,'zdufr zdvfr zdhfr'
494	c DO l=1,nlayer
495	c PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l)
496	c ENDDO
497	c STOP
498
499	DO l=1,nlayer
500	DO ig=1,ngrid
501	pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
502	pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
503	pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
504	ENDDO
505	ENDDO
506
507	DO ig=1,ngrid
508	zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig)
509	ENDDO
510
511	ELSE
512
513	DO j=jj_begin-offset,jj_end+offset
514	DO i=ii_begin-offset,ii_end+offset
515	ig=(j-1)*iim+i
516	zdtsrf(ig)=zdtsrf(ig)+
517	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
518	ENDDO
519	ENDDO
520
521	! write(66,*)"tsrf",maxval(zdtsrf(:)),minval(zdtsrf(:))
522	! write(66,*)"frd",maxval(fluxrad(:)),minval(fluxrad(:))
523	! write(66,*)"fgd",maxval(fluxgrd(:)),minval(fluxgrd(:))
524	ENDIF
525	c-----------------------------------------------------------------------
526	c 4. Dry convective adjustment:
527	c -----------------------------
528
529	IF(calladj) THEN
530
531	DO l=1,nlayer
532	DO ig=1,ngrid
533	zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l)
534	ENDDO
535	ENDDO
536	CALL zerophys(ngrid*nlayer,zdufr)
537	CALL zerophys(ngrid*nlayer,zdvfr)
538	CALL zerophys(ngrid*nlayer,zdhfr)
539	CALL convadj(ngrid,nlayer,ptimestep,
540	$ pplay,pplev,zpopsk,
541	$ pu,pv,zh,
542	$ pdu,pdv,zdum1,
543	$ zdufr,zdvfr,zdhfr)
544
545	DO l=1,nlayer
546	DO ig=1,ngrid
547	pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)
548	pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)
549	pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)
550	ENDDO
551	ENDDO
552
553	ENDIF
554
555	!101 continue
556	c-----------------------------------------------------------------------
557	c On incremente les tendances physiques de la temperature du sol:
558	c ---------------------------------------------------------------
559	! WRITE(55,*)"tsurf",maxval(tsurf(:)),minval(tsurf(:)),it
560
561	DO j=jj_begin-offset,jj_end+offset
562	DO i=ii_begin-offset,ii_end+offset
563	ig=(j-1)*iim+i
564	tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)
565	ENDDO
566	ENDDO
567	c-----------------------------------------------------------------------
568	c soil temperatures:
569	c --------------------
570
571	IF (callsoil) THEN
572	CALL soil(ngrid,nsoilmx,.false.,inertie,
573	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
574	IF(lverbose) THEN
575	PRINT*,'Surface Heat capacity,conduction Flux, Ts,
576	s dTs, dt'
577	PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout),
578	s zdtsrf(igout),ptimestep
579	ENDIF
580	ENDIF
581
582	c-----------------------------------------------------------------------
583	c sorties:
584	c --------
585	if(zday.GT.zday_last+period_sort) then
586	zday_last=zday
587	c Ecriture/extension de la coordonnee temps
588
589	do ig=1,ngrid
590	zpmer(ig)=pplev(ig,1)exp(pphi(ig,1)/(kappacpp*285.))
591	enddo
592	endif
593	c-----------------------------------------------------------------------
594	IF(lastcall) THEN
595	PRINT*,'Ecriture du fichier de reinitialiastion de la physique'
596	ENDIF
597
598	icount=icount+1
599	! RETURN
600	END SUBROUTINE phyparam_lmd
601	END MODULE PHY

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

Download in other formats: