Context Navigation

radiation_mod.F @ 149

Last change on this file since 149 was 149, checked in by sdubey, 11 years ago

Added few new routines to read NC files and compute diagnostics to r145.

Few routines of dry physics including radiation module, surface process and convective adjustment in new routine phyparam.f90. dynetat to read start files for dynamics.

check_conserve routine to compute conservation of quatities like mass, energy etc.etat0_heldsz.f90 for held-suarez test case initial conditions.

new Key time_style=lmd or dcmip to use day_step, ndays like in LMDZ

File size: 27.4 KB

Line
1	MODULE RADIATION
2	USE ICOSA
3	USE dimphys_mod
4	! USE PHY
5	contains
6
7	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8	SUBROUTINE zerophys(n,x)
9	IMPLICIT NONE
10	INTEGER n
11	REAL x(n)
12	INTEGER i
13
14	DO i=1,n
15	x(i)=0.
16	ENDDO
17	RETURN
18	END subroutine zerophys
19	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
20
21	SUBROUTINE solarlong(pday,psollong)
22	IMPLICIT NONE
23	c=======================================================================
24	c
25	c Objet:
26	c ------
27	c
28	c Calcul de la distance soleil-planete et de la declinaison
29	c en fonction du jour de l'annee.
30	c
31	c
32	c Methode:
33	c --------
34	c
35	c Calcul complet de l'elipse
36	c
37	c Interface:
38	c ----------
39	c
40	c Uncommon comprenant les parametres orbitaux.
41	c
42	c Arguments:
43	c ----------
44	c
45	c Input:
46	c ------
47	c pday jour de l'annee (le jour 0 correspondant a l'equinoxe)
48	c lwrite clef logique pour sorties de controle
49	c
50	c Output:
51	c -------
52	c pdist_sol distance entre le soleil et la planete
53	c ( en unite astronomique pour utiliser la constante
54	c solaire terrestre 1370 Wm-2 )
55	c pdecli declinaison ( en radians )
56	c
57	c=======================================================================
58	c-----------------------------------------------------------------------
59	c Declarations:
60	c -------------
61
62	c arguments:
63	c ----------
64
65	REAL pday,pdist_sol,pdecli,psollong
66	LOGICAL lwrite
67
68	c Local:
69	c ------
70
71	REAL zanom,xref,zx0,zdx,zteta,zz
72	INTEGER iter
73
74
75	c-----------------------------------------------------------------------
76	c calcul de l'angle polaire et de la distance au soleil :
77	c -------------------------------------------------------
78
79	c calcul de l'zanomalie moyenne
80
81	zz=(pday-peri_day)/year_day
82	zanom=2.pi(zz-nint(zz))
83	xref=abs(zanom)
84
85	c resolution de lequation horaire zx0 - e * sin (zx0) = xref
86	c methode de Newton
87
88	zx0=xref+e_elips*sin(xref)
89	DO 110 iter=1,10
90	zdx=-(zx0-e_elipssin(zx0)-xref)/(1.-e_elipscos(zx0))
91	if(abs(zdx).le.(1.e-7)) goto 120
92	zx0=zx0+zdx
93	110 continue
94	120 continue
95	zx0=zx0+zdx
96	if(zanom.lt.0.) zx0=-zx0
97
98	c zteta est la longitude solaire
99
100	zteta=2.atan(sqrt((1.+e_elips)/(1.-e_elips))tan(zx0/2.))
101
102	psollong=zteta-timeperi
103
104	IF(psollong.LT.0.) psollong=psollong+2.*pi
105	IF(psollong.GT.2.pi) psollong=psollong-2.pi
106
107	RETURN
108	END SUBROUTINE solarlong
109
110	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
111
112	SUBROUTINE orbite(pls,pdist_sol,pdecli)
113	IMPLICIT NONE
114	c====================================================================
115	c
116	c Objet:
117	c ------
118	c
119	c Distance from sun and declimation as a function of the solar
120	c longitude Ls
121	c
122	c Interface:
123	c ----------
124	c Arguments:
125	c ----------
126	c
127	c Input:
128	c ------
129	c pls Ls
130	c
131	c Output:
132	c -------
133	c pdist_sol Distance Sun-Planet in UA
134	c pdecli declinaison ( en radians )
135	c
136	c====================================================================
137	c-----------------------------------------------------------------------
138	c Declarations:
139	c -------------
140	c arguments:
141	c ----------
142
143	REAL pday,pdist_sol,pdecli,pls
144
145	c-----------------------------------------------------------------------
146
147	c Distance Sun-Planet
148
149	pdist_sol=p_elips/(1.+e_elips*cos(pls+timeperi))
150
151	c Solar declination
152
153	pdecli= asin (sin(pls)sin(obliquitpi/180.))
154
155	c-----------------------------------------------------------------------
156	c sorties eventuelles:
157	c ---------------------
158
159	END SUBROUTINE orbite
160	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
161
162	SUBROUTINE iniorbit
163	$ (paphelie,pperiheli,pyear_day,pperi_day,pobliq)
164	IMPLICIT NONE
165	c=====================================================================
166	c
167	c Auteur:
168	c -------
169	c Frederic Hourdin 22 Fevrier 1991
170	c
171	c Objet:
172	c ------
173	c Initialisation du sous programme orbite qui calcule
174	c a une date donnee de l'annee de duree year_day commencant
175	c a l'equinoxe de printemps et dont le perihelie se situe
176	c a la date peri_day, la distance au soleil et la declinaison.
177	c
178	c Interface:
179	c ----------
180	c - Doit etre appele avant d'utiliser orbite.
181	c - initialise le common comorbit
182	c
183	c Arguments:
184	c ----------
185	c
186	c Input:
187	c ------
188	c aphelie \ aphelie et perihelie de l'orbite
189	c periheli / en millions de kilometres.
190	c
191	c=====================================================================
192	c Declarations:
193	c -------------
194	c Arguments:
195	c ----------
196
197	REAL paphelie,pperiheli,pyear_day,pperi_day,pobliq
198
199	c Local:
200	c ------
201
202	REAL zxref,zanom,zz,zx0,zdx
203	INTEGER iter
204
205	c'-----------------------------------------------------------------------
206
207	aphelie =paphelie
208	periheli=pperiheli
209	year_day=pyear_day
210	obliquit=pobliq
211	peri_day=pperi_day
212
213	PRINT*,'Perihelie en Mkm ',periheli
214	PRINT*,'Aphelise en Mkm ',aphelie
215	PRINT*,'obliquite en degres :',obliquit
216	unitastr=149.597927
217	e_elips=(aphelie-periheli)/(periheli+aphelie)
218	p_elips=0.5(periheli+aphelie)(1-e_elips*e_elips)/unitastr
219
220	print*,'e_elips',e_elips
221	print*,'p_elips',p_elips
222
223	c-----------------------------------------------------------------------
224	c calcul de l'angle polaire et de la distance au soleil :
225	c -------------------------------------------------------
226
227	c calcul de l'zanomalie moyenne
228
229	zz=(year_day-pperi_day)/year_day
230	zanom=2.pi(zz-nint(zz))
231	zxref=abs(zanom)
232	PRINT*,'zanom ',zanom
233
234	c resolution de lequation horaire zx0 - e * sin (zx0) = zxref
235	c methode de Newton
236
237	zx0=zxref+e_elips*sin(zxref)
238	DO 110 iter=1,100
239	zdx=-(zx0-e_elipssin(zx0)-zxref)/(1.-e_elipscos(zx0))
240	if(abs(zdx).le.(1.e-12)) goto 120
241	zx0=zx0+zdx
242	110 continue
243	120 continue
244	zx0=zx0+zdx
245	if(zanom.lt.0.) zx0=-zx0
246	PRINT*,'zx0 ',zx0
247
248	c zteta est la longitude solaire
249
250	timeperi=2.atan(sqrt((1.+e_elips)/(1.-e_elips))tan(zx0/2.))
251	PRINT*,'longitude solaire du perihelie timeperi = ',timeperi
252
253	RETURN
254	END SUBROUTINE iniorbit
255	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
256
257	SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad)
258	IMPLICIT NONE
259
260	c=======================================================================
261	c
262	c Calcul of equivalent solar angle and and fraction of day whithout
263	c diurnal cycle.
264	c
265	c parmeters :
266	c -----------
267	c
268	c Input :
269	c -------
270	c npts number of points
271	c pdeclin solar declinaison
272	c plat(npts) latitude
273	c phaut hauteur typique de l'atmosphere
274	c prad rayon planetaire '
275	c
276	c Output :
277	c --------
278	c pmu(npts) equivalent cosinus of the solar angle
279	c pfract(npts) fractionnal day
280	c
281	c=======================================================================
282
283	c-----------------------------------------------------------------------
284	c
285	c 0. Declarations :
286	c -----------------
287
288	c Arguments :
289	c -----------
290	INTEGER npts
291	REAL plat(npts),pmu(npts), pfract(npts)
292	REAL phaut,prad,pdeclin
293	c
294	c Local variables :
295	c -----------------
296	INTEGER j,i,ij,ig
297	REAL pi
298	REAL z,cz,sz,tz,phi,cphi,sphi,tphi
299	REAL ap,a,t,b
300	REAL alph
301
302	c-----------------------------------------------------------------------
303
304	!print*,'npts,pdeclin'
305	!print*,npts,pdeclin
306	pi = 4. * atan(1.0)
307	!print*,'PI=',pi
308	pi=2.*asin(1.)
309	z = pdeclin
310	cz = cos (z)
311	sz = sin (z)
312	!print*,'cz,sz',cz,sz
313
314	! DO j=jj_begin-offset,jj_end+offset
315	! DO i=ii_begin-offset,ii_end+offset
316	! ig=(j-1)*iim+i
317
318	DO ig=1,npts
319	phi = plat(ig)
320	cphi = cos(phi)
321	if (cphi.le.1.e-9) cphi=1.e-9
322	sphi = sin(phi)
323	tphi = sphi / cphi
324	b = cphi * cz
325	t = -tphi * sz / cz
326	a = 1.0 - t*t
327	ap = a
328	IF(t.eq.0.) then
329	t=0.5*pi
330	ELSE
331	IF (a.lt.0.) a = 0.
332	t = sqrt(a) / t
333	IF (t.lt.0.) then
334	t = -atan (-t) + pi
335	ELSE
336	t = atan(t)
337	ENDIF
338	ENDIF
339
340	pmu(ig) = (sphiszt) / pi + b*sin(t)/pi
341	pfract(ig) = t / pi
342	IF (ap .lt.0.) then
343	pmu(ig) = sphi * sz
344	pfract(ig) = 1.0
345	ENDIF
346	IF (pmu(ig).le.0.0) pmu(ig) = 0.0
347	pmu(ig) = pmu(ig) / pfract(ig)
348	IF (pmu(ig).eq.0.) pfract(ig) = 0.
349	ENDDO
350	! ENDDO
351	c-----------------------------------------------------------------------
352	c correction de rotondite:
353	c ------------------------
354
355	alph=phaut/prad
356	! DO j=jj_begin-offset,jj_end+offset
357	! DO i=ii_begin-offset,ii_end+offset
358	! ig=(j-1)*iim+i
359	DO ig = 1,npts
360	pmu(ig)=sqrt(1224.pmu(ig)pmu(ig)+1.)/35.
361
362	c $ (sqrt(alphalphpmu(ig)pmu(ig)+2.alph+1.)-alph*pmu(ig))
363	ENDDO
364	! ENDDO
365
366	RETURN
367	END SUBROUTINE mucorr
368	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
369
370	SUBROUTINE sw(ngrid,nlayer,ldiurn,
371	$ coefvis,albedo,
372	$ plevel,ps_rad,pmu,pfract,psolarf0,
373	$ fsrfvis,dtsw,
374	$ lwrite)
375	IMPLICIT NONE
376	c=======================================================================
377	c
378	c Rayonnement solaire en atmosphere non diffusante avec un
379	c coefficient d'absoprption gris.
380	c'
381	c=======================================================================
382	c
383	c declarations:
384	c -------------
385	c arguments:
386	c ----------
387	c
388	INTEGER ngrid,nlayer,i,j,ij
389	REAL albedo(ngrid),coefvis
390	REAL pmu(ngrid),pfract(ngrid)
391	REAL plevel(ngrid,nlayer+1),ps_rad
392	REAL psolarf0
393	REAL fsrfvis(ngrid),dtsw(ngrid,nlayer)
394	LOGICAL lwrite,ldiurn
395	c
396	c variables locales:
397	c ------------------
398	c
399
400	REAL zalb(ngrid),zmu(ngrid),zfract(ngrid)
401	REAL zplev(ngrid,nlayer+1)
402	REAL zflux(ngrid),zdtsw(ngrid,nlayer)
403
404	INTEGER ig,l,nlevel,ncount,igout
405	INTEGER,DIMENSION(ngrid)::index
406	REAL ztrdir(ngrid,nlayer+1),ztrref(ngrid,nlayer+1)
407	REAL zfsrfref(ngrid)
408	REAL z1(ngrid)
409	REAL zu(ngrid,nlayer+1)
410	REAL tau0
411
412	EXTERNAL SSUM
413	EXTERNAL ismax
414	REAL ismax
415
416	LOGICAL firstcall
417	SAVE firstcall
418	DATA firstcall/.true./
419
420	c-----------------------------------------------------------------------
421	c 1. initialisations:
422	c -------------------
423
424
425
426	IF (firstcall) THEN
427
428	IF (ngrid.NE.ngrid) THEN
429	PRINT*,'STOP in inifis'
430	PRINT*,'Probleme de dimenesions :'
431	PRINT*,'ngrid = ',ngrid
432	PRINT*,'ngrid = ',ngrid
433	STOP
434	ENDIF
435
436
437	IF (nlayer.NE.llm) THEN
438	PRINT*,'STOP in inifis'
439	PRINT*,'Probleme de dimenesions :'
440	PRINT*,'nlayer = ',nlayer
441	PRINT*,'llm = ',llm
442	STOP
443	ENDIF
444
445	ENDIF
446
447	nlevel=nlayer+1
448	c-----------------------------------------------------------------------
449	c Definitions des tableaux locaux pour les points ensoleilles:
450	c ------------------------------------------------------------
451
452	IF (ldiurn) THEN
453	ncount=0
454	DO j=jj_begin-offset,jj_end+offset
455	DO i=ii_begin-offset,ii_end+offset
456	ig=(j-1)*iim+i
457	index(ig)=0
458	ENDDO
459	ENDDO
460	DO j=jj_begin-offset,jj_end+offset
461	DO i=ii_begin-offset,ii_end+offset
462	ig=(j-1)*iim+i
463	IF(pfract(ig).GT.1.e-6) THEN
464	ncount=ncount+1
465	index(ncount)=ig
466	ENDIF
467	ENDDO
468	ENDDO
469	! SARVESH CHANGED NCOUNT TO NGRID TO BE CONSISTENT ???
470
471	CALL monGATHER(ncount,zfract,pfract,index)
472	CALL monGATHER(ncount,zmu,pmu,index)
473	CALL monGATHER(ncount,zalb,albedo,index)
474	DO l=1,nlevel
475	CALL monGATHER(ncount,zplev(1,l),plevel(1,l),index)
476	ENDDO
477	ELSE
478	ncount=ngrid
479	zfract(:)=pfract(:)
480	zmu(:)=pmu(:)
481	zalb(:)=albedo(:)
482	zplev(:,:)=plevel(:,:)
483	ENDIF
484
485	c-----------------------------------------------------------------------
486	c calcul des profondeurs optiques integres depuis p=0:
487	c ----------------------------------------------------
488
489	tau0=-.5*log(coefvis)
490
491	c calcul de la partie homogene de l'opacite
492	c'
493	tau0=tau0/ps_rad
494
495
496	DO l=1,nlayer+1
497	DO j=jj_begin-offset,jj_end+offset
498	DO i=ii_begin-offset,ii_end+offset
499	ig=(j-1)*iim+i
500	zu(ig,l)=tau0*zplev(ig,l)
501	ENDDO
502	ENDDO
503	ENDDO
504
505	c-----------------------------------------------------------------------
506	c 2. calcul de la transmission depuis le sommet de l'atmosphere:
507	c' -----------------------------------------------------------
508
509	DO l=1,nlevel
510	DO j=jj_begin-offset,jj_end+offset
511	DO i=ii_begin-offset,ii_end+offset
512	ig=(j-1)*iim+i
513	ztrdir(ig,l)=exp(-zu(ig,l)/zmu(ig))
514	ENDDO
515	ENDDO
516	ENDDO
517
518	IF (lwrite) THEN
519	igout=ncount/2+1
520	PRINT*
521	PRINT*,'Diagnostique des transmission dans le spectre solaire'
522	PRINT*,'zfract, zmu, zalb'
523	PRINT*,zfract(igout), zmu(igout), zalb(igout)
524	PRINT*,'Pression, quantite d abs, transmission'
525	DO l=1,nlayer+1
526	PRINT*,zplev(igout,l),zu(igout,l),ztrdir(igout,l)
527	ENDDO
528	ENDIF
529
530	c-----------------------------------------------------------------------
531	c 3. taux de chauffage, ray. solaire direct:
532	c ------------------------------------------
533
534	DO l=1,nlayer
535	DO j=jj_begin-offset,jj_end+offset
536	DO i=ii_begin-offset,ii_end+offset
537	ig=(j-1)*iim+i
538	zdtsw(ig,l)=gpsolarf0zfract(ig)zmu(ig)
539	$ (ztrdir(ig,l+1)-ztrdir(ig,l))/
540	$ (cpp*(zplev(ig,l)-zplev(ig,l+1)))
541	ENDDO
542	ENDDO
543	ENDDO
544	IF (lwrite) THEN
545	PRINT*
546	PRINT*,'Diagnostique des taux de chauffage solaires:'
547	PRINT*,' 1 taux de chauffage lie au ray. solaire direct'
548	DO l=1,nlayer
549	PRINT*,zdtsw(igout,l)
550	ENDDO
551	ENDIF
552
553
554	c-----------------------------------------------------------------------
555	c 4. calcul du flux solaire arrivant sur le sol:
556	c ----------------------------------------------
557
558	DO j=jj_begin-offset,jj_end+offset
559	DO i=ii_begin-offset,ii_end+offset
560	ig=(j-1)*iim+i
561	z1(ig)=zfract(ig)zmu(ig)psolarf0*ztrdir(ig,1)
562	zflux(ig)=(1.-zalb(ig))*z1(ig)
563	zfsrfref(ig)= zalb(ig)*z1(ig)
564	ENDDO
565	ENDDO
566	IF (lwrite) THEN
567	PRINT*
568	PRINT*,'Diagnostique des taux de chauffage solaires:'
569	PRINT*,' 2 flux solaire net incident sur le sol'
570	PRINT*,zflux(igout)
571	ENDIF
572
573
574	c-----------------------------------------------------------------------
575	c 5.calcul des traansmissions depuis le sol, cas diffus:
576	c ------------------------------------------------------
577
578	DO l=1,nlevel
579	DO j=jj_begin-offset,jj_end+offset
580	DO i=ii_begin-offset,ii_end+offset
581	ig=(j-1)*iim+i
582	ztrref(ig,l)=exp(-(zu(ig,1)-zu(ig,l))*1.66)
583	ENDDO
584	ENDDO
585	ENDDO
586
587	IF (lwrite) THEN
588	PRINT*
589	PRINT*,'Diagnostique des taux de chauffage solaires'
590	PRINT*,' 3 transmission avec les sol'
591	PRINT*,'niveau transmission'
592	DO l=1,nlevel
593	PRINT*,l,ztrref(igout,l)
594	ENDDO
595	ENDIF
596
597	c-----------------------------------------------------------------------
598	c 6.ajout a l'echauffement de la contribution du ray. sol. reflechit:
599	c' -------------------------------------------------------------------
600
601	DO l=1,nlayer
602	DO j=jj_begin-offset,jj_end+offset
603	DO i=ii_begin-offset,ii_end+offset
604	ig=(j-1)*iim+i
605	zdtsw(ig,l)=zdtsw(ig,l)+
606	$ gzfsrfref(ig)(ztrref(ig,l+1)-ztrref(ig,l))/
607	$ (cpp*(zplev(ig,l+1)-zplev(ig,l)))
608	ENDDO
609	ENDDO
610	ENDDO
611
612	c-----------------------------------------------------------------------
613	c 10. sorites eventuelles:
614	c ------------------------
615
616	IF (lwrite) THEN
617	PRINT*
618	PRINT*,'Diagnostique des taux de chauffage solaires:'
619	PRINT*,' 3 taux de chauffage total'
620	DO l=1,nlayer
621	PRINT*,zdtsw(igout,l)
622	ENDDO
623	ENDIF
624
625	IF (ldiurn) THEN
626	CALL zerophys(ngrid,fsrfvis)
627	CALL monscatter(ncount,fsrfvis,index,zflux)
628	CALL zerophys(ngrid*nlayer,dtsw)
629	DO l=1,nlayer
630	CALL monscatter(ncount,dtsw(1,l),index,zdtsw(1,l))
631	ENDDO
632	ELSE
633	!print*,'NOT DIURNE'
634	fsrfvis(:)=zflux(:)
635	dtsw(:,:)=zdtsw(:,:)
636	ENDIF
637
638	RETURN
639	END SUBROUTINE sw
640	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
641	SUBROUTINE lw(ngrid,nlayer,coefir,emissiv,
642	$ pp,ps_rad,ptsurf,pt,
643	$ pfluxir,pdtlw,
644	$ lwrite)
645
646	IMPLICIT NONE
647	c=======================================================================
648	c
649	c calcul de l'evolution de la temperature sous l'effet du rayonnement
650	c infra-rouge.
651	c Pour simplifier, les transmissions sont precalculees et ne
652	c dependent que de l'altitude.
653	c
654	c arguments:
655	c ----------
656	c'
657	c entree:
658	c -------
659	c ngrid nombres de points de la grille horizontale
660	c nlayer nombre de couches
661	c ptsurf(ngrid) temperature de la surface
662	c pt(ngrid,nlayer) temperature des couches
663	c pp(ngrid,nlayer+1) pression entre les couches
664	c lwrite variable logique pour sorties
665	c
666	c sortie:
667	c -------
668	c pdtlw(ngrid,nlayer) taux de refroidissement
669	c pfluxir(ngrid) flux infrarouge sur le sol
670	c
671	c=======================================================================
672
673	!c declarations:
674	!c -------------
675	!c arguments:
676	!c' ----------
677
678	INTEGER ngrid,nlayer
679	REAL coefir,emissiv(ngrid),ps_rad
680	REAL ptsurf(ngrid),pt(ngrid,nlayer),pp(ngrid,nlayer+1)
681	REAL pdtlw(ngrid,nlayer),pfluxir(ngrid)
682	LOGICAL lwrite
683
684	c variables locales:
685	c ------------------
686
687	INTEGER nlevel,ilev,ig,i,il
688	REAL zplanck(ngridmx,nlayermx+1),zcoef
689	REAL zfluxup(ngridmx,nlayermx+1),zfluxdn(ngridmx,nlayermx+1)
690	REAL zflux(ngridmx,nlayermx+1)
691	REAL zlwtr1(ngridmx),zlwtr2(ngridmx)
692	REAL zup(ngridmx,nlayermx+1),zdup(ngridmx)
693	REAL stephan
694
695	LOGICAL lstrong
696	SAVE lstrong,stephan
697	DATA lstrong/.true./
698	c-----------------------------------------------------------------------
699	c initialisations:
700	c ----------------
701
702	nlevel=nlayer+1
703	stephan=5.67e-08
704
705
706	pfluxir=0.0
707	pdtlw=0.0
708	!print*,"ngr,nlay,coefi",ngrid,nlayer,coefir
709	c-----------------------------------------------------------------------
710	c 2. calcul des quantites d'absorbants:
711	c' -------------------------------------
712
713	c absorption forte
714	IF(lstrong) THEN
715	DO ilev=1,nlevel
716	DO ig=1,ngrid
717	zup(ig,ilev)=pp(ig,ilev)pp(ig,ilev)/(2.g)
718	ENDDO
719	ENDDO
720	IF(lwrite) THEN
721	DO ilev=1,nlayer
722	PRINT*,' up(',ilev,') = ',zup(ngrid/2+1,ilev)
723	ENDDO
724	ENDIF
725	zcoef=-log(coefir)/sqrt(ps_radps_rad/(2.g))
726
727	c absorption faible
728	ELSE
729	DO ilev=1,nlevel
730	DO ig=1,ngrid
731	zup(ig,ilev)=pp(ig,ilev)
732	ENDDO
733	ENDDO
734	zcoef=-log(coefir)/ps_rad
735	ENDIF
736
737
738	c-----------------------------------------------------------------------
739	c 2. calcul de la fonction de corps noir:
740	c ---------------------------------------
741
742	DO ilev=1,nlayer
743	DO ig=1,ngrid
744	zplanck(ig,ilev)=pt(ig,ilev)*pt(ig,ilev)
745	zplanck(ig,ilev)=stephan*
746	$ zplanck(ig,ilev)*zplanck(ig,ilev)
747	ENDDO
748	ENDDO
749
750	c-----------------------------------------------------------------------
751	c 4. flux descendants:
752	c --------------------
753
754	DO ilev=1,nlayer
755	DO ig=1,ngrid
756	zfluxdn(ig,ilev)=0.
757	ENDDO
758	DO ig=1,ngrid
759	zdup(ig)=zup(ig,ilev)-zup(ig,nlevel)
760	ENDDO
761	CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
762
763	DO il=nlayer,ilev,-1
764	zlwtr2(:)=zlwtr1(:)
765	DO ig=1,ngrid
766	zdup(ig)=zup(ig,ilev)-zup(ig,il)
767	ENDDO
768	CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
769	DO ig=1,ngrid
770	zfluxdn(ig,ilev)=zfluxdn(ig,ilev)+
771	$ zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
772	ENDDO
773	ENDDO
774	ENDDO
775
776	DO ig=1,ngrid
777	zfluxdn(ig,nlevel)=0.
778	pfluxir(ig)=emissiv(ig)*zfluxdn(ig,1)
779	ENDDO
780
781	DO ig=1,ngrid
782	zfluxup(ig,1)=ptsurf(ig)*ptsurf(ig)
783	zfluxup(ig,1)=emissiv(ig)stephanzfluxup(ig,1)*zfluxup(ig,1)
784	$ +(1.-emissiv(ig))*zfluxdn(ig,1)
785	ENDDO
786
787	c-----------------------------------------------------------------------
788	c 3. flux montants:
789	c ------------------
790
791	DO ilev=1,nlayer
792	DO ig=1,ngrid
793	zdup(ig)=zup(ig,1)-zup(ig,ilev+1)
794	ENDDO
795	CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
796	DO ig=1,ngrid
797	zfluxup(ig,ilev+1)=zfluxup(ig,1)*zlwtr1(ig)
798	ENDDO
799	DO il=1,ilev
800	zlwtr2(:)=zlwtr1(:)
801	DO ig=1,ngrid
802	zdup(ig)=zup(ig,il+1)-zup(ig,ilev+1)
803	ENDDO
804	CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
805	DO ig=1,ngrid
806	zfluxup(ig,ilev+1)=zfluxup(ig,ilev+1)+
807	$ zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
808	ENDDO
809	ENDDO
810
811	ENDDO
812
813	c-----------------------------------------------------------------------
814	c 5. calcul des flux nets:
815	c ------------------------
816
817	DO ilev=1,nlevel
818	DO ig=1,ngrid
819	zflux(ig,ilev)=zfluxup(ig,ilev)-zfluxdn(ig,ilev)
820	ENDDO
821	ENDDO
822
823	c-----------------------------------------------------------------------
824	c 6. Calcul des taux de refroidissement:
825	c --------------------------------------
826
827	DO ilev=1,nlayer
828	DO ig=1,ngrid
829	pdtlw(ig,ilev)=(zflux(ig,ilev+1)-zflux(ig,ilev))*
830	$ g/(cpp*(pp(ig,ilev+1)-pp(ig,ilev)))
831	ENDDO
832	ENDDO
833
834	c-----------------------------------------------------------------------
835	c 10. sorties eventuelles:
836	c ------------------------
837
838	IF (lwrite) THEN
839
840	PRINT*
841	PRINT*,'Diagnostique rayonnement thermique'
842	PRINT*
843	PRINT*,'temperature ',
844	$ 'flux montant flux desc. taux de refroid.'
845	i=ngrid/2+1
846	WRITE(6,9000) ptsurf(i)
847	DO ilev=1,nlayer
848	WRITE(6,'(i4,4e18.4)') ilev,pt(i,ilev),
849	$ zfluxup(i,ilev),zfluxdn(i,ilev),pdtlw(i,ilev)
850	ENDDO
851	WRITE(6,9000) zfluxup(i,nlevel),zfluxdn(i,nlevel)
852
853	ENDIF
854
855	c-----------------------------------------------------------------------
856
857	RETURN
858	9000 FORMAT(4e18.4)
859	END SUBROUTINE lw
860	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
861
862	subroutine solang ( kgrid,psilon,pcolon,psilat,pcolat,
863	& ptim1,ptim2,ptim3,pmu0,pfract )
864	IMPLICIT NONE
865
866	C
867	C**** LW - ORGANIZES THE LONGWAVE CALCULATIONS
868	C
869	C PURPOSE.
870	C --------
871	C CALCULATES THE SOLAR ANGLE FOR ALL THE POINTS OF THE GRID
872	C ==== INPUTS ===
873	C
874	C PSILON(KGRID) : SINUS OF THE LONGITUDE
875	C PCOLON(KGRID) : COSINUS OF THE LONGITUDE
876	C PSILAT(KGRID) : SINUS OF THE LATITUDE
877	C PCOLAT(KGRID) : COSINUS OF THE LATITUDE
878	C PTIM1 : SIN(DECLI)
879	C PTIM2 : COS(DECLI)*COS(TIME)
880	C PTIM3 : SIN(DECLI)*SIN(TIME)
881	C
882	C ==== OUTPUTS ===
883	C
884	C PMU0 (KGRID) : SOLAR ANGLE
885	C PFRACT(KGRID) : DAY FRACTION OF THE TIME INTERVAL
886	C
887	C IMPLICIT ARGUMENTS : NONE
888	C --------------------
889	C
890	C METHOD.
891	C -------
892	C
893	C EXTERNALS.
894	C ----------
895	C
896	C NONE
897	C
898	C REFERENCE.
899	C ----------
900	C
901	C RADIATIVE PROCESSES IN METEOROLOGIE AND CLIMATOLOGIE
902	C PALTRIDGE AND PLATT
903	C
904	C AUTHOR.
905	C -------
906	C FREDERIC HOURDIN
907	C
908	C MODIFICATIONS.
909	C --------------
910	C ORIGINAL :90-01-14
911	C 92-02-14 CALCULATIONS DONE THE ENTIER GRID (J.Polcher)
912	C-----------------------------------------------------------------------
913	C
914	C ------------------------------------------------------------------
915
916	C-----------------------------------------------------------------------
917	C
918	C* 0.1 ARGUMENTS
919	C ---------
920	C
921	INTEGER kgrid
922	REAL ptim1,ptim2,ptim3
923	REAL psilon(kgrid),pcolon(kgrid),pmu0(kgrid),pfract(kgrid)
924	REAL psilat(kgrid), pcolat(kgrid)
925	C
926	INTEGER jl,i,j
927	REAL ztim1,ztim2,ztim3
928	C------------------------------------------------------------------------
929	C------------------------------------------------------------------------
930	C---------------------------------------------------------------------
931	C
932	C--------------------------------------------------------------------
933	C
934	C* 1. INITIALISATION
935	C --------------
936	C
937	!!!!!! 100 CONTINUE
938	C
939	DO j=jj_begin-offset,jj_end+offset
940	DO i=ii_begin-offset,ii_end+offset
941	jl=(j-1)*iim+i
942	pmu0(jl)=0.
943	pfract(jl)=0.
944	ENDDO
945	ENDDO
946	!C
947	!C* 1.1 COMPUTATION OF THE SOLAR ANGLE
948	!C ------------------------------
949	!C
950	DO j=jj_begin-offset,jj_end+offset
951	DO i=ii_begin-offset,ii_end+offset
952	jl=(j-1)*iim+i
953	ztim1=psilat(jl)*ptim1
954	ztim2=pcolat(jl)*ptim2
955	ztim3=pcolat(jl)*ptim3
956	pmu0(jl)=ztim1+ztim2pcolon(jl)+ztim3psilon(jl)
957	ENDDO
958	ENDDO
959	!C
960	!C* 1.2 DISTINCTION BETWEEN DAY AND NIGHT
961	!C ---------------------------------
962	!C
963	DO j=jj_begin-offset,jj_end+offset
964	DO i=ii_begin-offset,ii_end+offset
965	jl=(j-1)*iim+i
966	IF (pmu0(jl).gt.0.) THEN
967	pfract(jl)=1.
968	ELSE
969	pmu0(jl)=0.
970	pfract(jl)=0.
971	ENDIF
972	ENDDO
973	ENDDO
974	RETURN
975	END SUBROUTINE solang
976	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
977
978	SUBROUTINE monGATHER(n,a,b,index)
979	c
980	IMPLICIT NONE
981	C
982	INTEGER n,ng,index(n),i,j,ij
983	REAL a(n),b(n)
984	c
985	DO 100 i=1,n
986	a(i)=b(index(i))
987	100 CONTINUE
988
989	! DO j=jj_begin-offset,jj_end+offset
990	! DO i=ii_begin-offset,ii_end+offset
991	! ij=(j-1)*iim+i
992	! a(ij)=b(index(ij))
993	!c
994	RETURN
995	END SUBROUTINE monGATHER
996	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
997
998	subroutine monscatter(n,a,index,b)
999	C
1000	implicit none
1001	integer N,INDEX(n),I
1002	real A(n),B(n)
1003	c
1004	c
1005	DO 100 I=1,N
1006	A(INDEX(I))=B(I)
1007	100 CONTINUE
1008	c
1009	return
1010	end SUBROUTINE monscatter
1011	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1012
1013	SUBROUTINE lwtr(ngrid,coef,lstrong,dup,transm)
1014	IMPLICIT NONE
1015	INTEGER ngrid
1016	REAL coef
1017	LOGICAL lstrong
1018	REAL dup(ngrid),transm(ngrid)
1019	INTEGER ig
1020
1021	IF(lstrong) THEN
1022	DO ig=1,ngrid
1023	transm(ig)=exp(-coef*sqrt(dup(ig)))
1024	ENDDO
1025	ELSE
1026	DO ig=1,ngrid
1027	transm(ig)=exp(-coef*dup(ig))
1028	ENDDO
1029	ENDIF
1030	RETURN
1031	END subroutine lwtr
1032	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1033	END MODULE RADIATION

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

Context Navigation

source: codes/icosagcm/trunk/src/radiation_mod.F @ 149

Download in other formats: