/[MITgcm]/MITgcm/pkg/aim/phy_radiat.F
ViewVC logotype

Diff of /MITgcm/pkg/aim/phy_radiat.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by cnh, Fri Jan 26 00:14:32 2001 UTC revision 1.2 by adcroft, Fri Feb 2 21:36:29 2001 UTC
# Line 0  Line 1 
1    C $Header$
2    C $Name$
3    
4          SUBROUTINE SOL_OZ (SOLC,TYEAR,FSOL,OZONE)
5    
6    C--
7    C--   SUBROUTINE SOL_OZ (SOLC,TYEAR,FSOL,OZONE)
8    C--
9    C--   Purpose: Compute the flux of incoming solar radiation
10    C--            and a climatological ozone profile for SW absorption
11    C--   Input:   SOLC   = solar constant (area averaged)
12    C--            TYEAR  = time as fraction of year (0-1, 0 = 1jan.h00)
13    C--   Output:  FSOL   = flux of incoming solar radiation        (2-dim)
14    C--            OZONE  = strat. ozone as fraction of global mean (2-dim)
15    C--
16    
17    
18          IMPLICIT rEAL*8 (A-H,O-Z)
19    
20    
21    C     Resolution parameters
22    C
23    #include "atparam.h"
24    #include "atparam1.h"
25    C
26          PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT )
27    C
28    C     Constants + functions of sigma and latitude
29    C
30    #include "com_physcon.h"
31    
32          REAL FSOL(NLON,NLAT), OZONE(NLON,NLAT)
33    C
34    C     ALPHA = year phase ( 0 - 2pi, 0 = winter solstice = 22dec.h00 )
35          ALPHA=4.*ASIN(1.)*(TYEAR+10./365.)
36    
37          CSR1=-0.796*COS(ALPHA)
38          CSR2= 0.147*COS(2*ALPHA)-0.477
39          COZ1= 0.0
40    C     COZ1= 0.2*SIN(ALPHA)
41          COZ2= 0.3
42    C
43          DO J=1,NLAT
44            FSOL(1,J)=SOLC*MAX(0.,1.0+CSR1*FMU(J,1)+CSR2*FMU(J,2))
45            OZONE(1,J)=1.0+COZ1*FMU(J,1)+COZ2*FMU(J,2)
46            DO I=2,NLON
47              FSOL(I,J)=FSOL(1,J)
48              OZONE(I,J)=OZONE(1,J)
49            ENDDO
50          ENDDO
51    C
52          RETURN
53          END
54    
55    
56          SUBROUTINE RADSW (PSA,QA,RH,
57         *                  FSOL,OZONE,ALB,
58         *                  CLOUDC,FTOP,FSFC,DFABS)
59    C--
60    C--  SUBROUTINE RADSW (PSA,QA,RH,
61    C-- *                  FSOL,OZONE,ALB,
62    C-- *                  CLOUDC,FTOP,FSFC,DFABS)
63    C--
64    C--   Purpose: Compute the absorption of shortwave radiation and
65    C--            initialize arrays for longwave-radiation routines
66    C--   Input:   PSA    = norm. surface pressure [p/p0]           (2-dim)
67    C--            QA     = specific humidity [g/kg]                (3-dim)
68    C--            RH     = relative humidity                       (3-dim)
69    C--            FSOL   = flux of incoming solar radiation        (2-dim)
70    C--            OZONE  = strat. ozone as fraction of global mean (2-dim)
71    C--            ALB    = surface albedo                          (2-dim)
72    C--   Output:  CLOUDC = total cloud cover                           (2-dim)
73    C--            FTOP   = net downw. flux of sw rad. at the atm. top  (2-dim)
74    C--            FSFC   = net downw. flux of sw rad. at the surface   (2-dim)
75    C--            DFABS  = flux of sw rad. absorbed by each atm. layer (3-dim)
76    C--
77    
78    
79          IMPLICIT rEAL*8 (A-H,O-Z)
80    
81    
82    C     Resolution parameters
83    C
84    #include "atparam.h"
85    #include "atparam1.h"
86    #include "Lev_def.h"
87    C
88          PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT )
89    C
90    C     Constants + functions of sigma and latitude
91    C
92    #include "com_physcon.h"
93    C
94    C     Radiation parameters
95    C
96    #include "com_radcon.h"
97    C
98          REAL PSA(NGP), QA(NGP,NLEV), RH(NGP,NLEV),
99         *     FSOL(NGP), OZONE(NGP), ALB(NGP)
100    
101          REAL CLOUDC(NGP), FTOP(NGP), FSFC(NGP), DFABS(NGP,NLEV)
102    
103          REAL FLUX(NGP), FREFL(NGP), TAUOZ(NGP)
104          INTEGER NL1(NGP)
105    Cchdbg
106          INTEGER Npas
107          SAVE npas
108          LOGICAL Ifirst
109          SAVE Ifirst
110          DATA Ifirst /.TRUE./
111          REAL clsum(NGP)
112          SAVE clsum
113          REAL ABWLW1
114    cchdbg
115    C
116          DO J=1,NGP
117            NL1(J)=NLEVxy(J)-1
118          ENDDO
119    C
120    C--   1.  Cloud cover:
121    C         defined as a linear fun. of the maximum relative humidity
122    C         in all tropospheric layers above PBL:
123    C         CLOUDC =  0 for RHmax < RHCL1, = 1 for RHmax > RHCL2.
124    C         This value is reduced by a factor (Qbase/QACL) if the
125    C         cloud-base absolute humidity Qbase < QACL.
126    C
127          DRHCL=RHCL2-RHCL1
128          RCL=1./(DRHCL*QACL)
129    C
130          DO 122 J=1,NGP
131            CLOUDC(J)=0.
132      122 CONTINUE
133    C
134          DO 123 K=1,NLEV
135          DO 123 J=1,NGP
136            DFABS(J,K)=0.
137      123 CONTINUE
138    C
139          DO 124 J=1,NGP
140          DO 124 K=2,NL1(J)
141            CLOUDC(J)=MAX(CLOUDC(J),(RH(J,K)-RHCL1))
142      124 CONTINUE
143    C
144          DO 126 J=1,NGP
145            IF ( NL1(J) .GT. 0 ) THEN
146             CLOUDC(J)=MIN(CLOUDC(J),DRHCL)*MIN(QA(J,NL1(J)),QACL)*RCL
147            ENDIF
148    cchdbg *******************************************
149    cchdbg        CLOUDC(J)=MIN(CLOUDC(J),DRHCL)/DRHCL
150    cchdbg *******************************************
151    clear sky experiment
152    c       cloudc(j) = 0.
153      126 CONTINUE
154    C
155    C
156    C--   2. Shortwave transmissivity:
157    C        function of layer mass, ozone (in the statosphere),
158    C        abs. humidity and cloud cover (in the troposphere)
159    C
160          DO 202 J=1,NGP
161            TAU(J,1)=EXP(-ABSSW*PSA(J)*DSIG(1))
162            TAUOZ(J)=EXP(-EPSSW*OZONE(J)*PSA(J))
163      202 CONTINUE
164    C
165    chhh      WRITE(0,*) ' Hello from RADSW'
166          DO 204 J=1,NGP
167          DO 204 K=2,NL1(J)
168            TAU(J,K)=EXP(-(ABSSW+ABWSW*QA(J,K)
169         *           +ABCSW*CLOUDC(J)*QA(J,NL1(J)))*PSA(J)*DSIG(K))
170      204 CONTINUE
171    
172          DO 206 J=1,NGP
173           IF ( NLEVxy(J) .GT. 0 ) THEN
174            TAU(J,NLEVxy(J))=EXP(-(ABSSW+ABWSW*QA(J,NLEVxy(J)))
175         *                              *PSA(J)*DSIG(NLEVxy(J)))
176           ENDIF
177      206 CONTINUE
178    C
179    C---  3. Shortwave downward flux
180    C      
181    C     3.1  Absorption in the stratosphere
182    C
183          DO 312 J=1,NGP
184            FLUX(J)=TAU(J,1)*TAUOZ(J)*FSOL(J)
185            DFABS(J,1)=FSOL(J)-FLUX(J)
186      312 CONTINUE
187    C      
188    C     3.2  Reflection at the top of the troposphere
189    C          (proportional to cloud cover).
190    C
191          DO 322 J=1,NGP
192            FREFL(J)=ALBCL*CLOUDC(J)*FLUX(J)
193            FTOP(J) =FSOL(J)-FREFL(J)
194            FLUX(J) =FLUX(J)-FREFL(J)
195      322 CONTINUE  
196    C      
197    C     3.3  Absorption in the troposphere
198    C
199          DO 332 J=1,NGP
200          DO 332 K=2,NLEVxy(J)
201            DFABS(J,K)=FLUX(J)
202            FLUX(J)=TAU(J,K)*FLUX(J)
203            DFABS(J,K)=DFABS(J,K)-FLUX(J)
204      332 CONTINUE
205    C
206    C---  4. Shortwave upward flux
207    C      
208    C     4.1  Absorption and reflection at the surface
209    C
210          DO 412 J=1,NGP
211            FREFL(J)=ALB(J)*FLUX(J)
212            FSFC(J) =FLUX(J)-FREFL(J)
213            FLUX(J) =FREFL(J)
214      412 CONTINUE
215    C      
216    C     4.2  Absorption in the atmosphere
217    C
218          DO 422 J=1,NGP
219          DO 422 K=NLEVxy(J),1,-1
220            DFABS(J,K)=DFABS(J,K)+FLUX(J)
221            FLUX(J)=TAU(J,K)*FLUX(J)
222            DFABS(J,K)=DFABS(J,K)-FLUX(J)
223      422 CONTINUE
224    C
225    C     4.3  Absorbed solar radiation = incoming - outgoing
226    C
227          DO 432 J=1,NGP
228            FTOP(J)=FTOP(J)-FLUX(J)
229      432 CONTINUE
230    cdj
231    c     write(0,*)'position j=20'
232    c     j=20
233    c     write(0,*)'ftop fsfc ftop-fsfc'
234    c     write(0,*)ftop(j),fsfc(j),ftop(j)-fsfc(j)
235    c     write(0,*)
236    c     write(0,*)'k dfabs'
237    c     do k = 1, nlevxy(j)
238    c       write(0,*)k,dfabs(j,k)
239    c     enddo
240    c     write(0,*)'sum dfabs'
241    c     write(0,*)sum(dfabs(j,:))
242    cdj
243    C
244    C---  5.  Initialization of longwave radiation model
245    C
246    C     5.1  Longwave transmissivity:
247    C          function of layer mass, abs. humidity and cloud cover.
248    C
249          DO 512 J=1,NGP
250            TAU(J,1)=EXP(-ABSLW*PSA(J)*DSIG(1))
251      512 CONTINUE
252    C
253          DO 514 J=1,NGP
254          DO 514 K=2,NL1(J)
255            TAU(J,K)=EXP(-(ABSLW+ABWLW*QA(J,K)
256         *           +ABCLW*CLOUDC(J)*QA(J,NL1(J)))*PSA(J)*DSIG(K))
257      514 CONTINUE
258    C
259    cchdbg ***************************************************
260    c      ABCLW1=0.15
261    c      DO 514 J=1,NGP
262    c      DO 514 K=2,NL1(J)-1
263    c        TAU(J,K)=EXP(-(ABSLW+ABWLW*QA(J,K)
264    c     *           +ABCLW1*CLOUDC(J)*QA(J,NL1(J)))*PSA(J)*DSIG(K))
265    c  514 CONTINUE
266    C
267    c      DO 515 J=1,NGP
268    c      DO 515 K=NL1(J),NL1(J)
269    c        TAU(J,K)=EXP(-(ABSLW+ABWLW*QA(J,K)
270    c     *           +ABCLW*CLOUDC(J))*PSA(J)*DSIG(K))
271    c  515 CONTINUE
272    cchdbg ************************************************************
273    C
274    C *********************************************************************
275    C *********************************************************************
276    C *****************************************************************
277    cchdbg
278    c      if(Ifirst) then
279    c        npas=0
280    c        do J=1,NGP
281    c          clsum(J)=0.
282    c        enddo
283    c        ifirst=.FALSE.
284    c      ENDIF
285    C
286    c      npas=npas+1
287    c      DO J=1,NGP
288    c        clsum(J)=clsum(J)+ABCLW*CLOUDC(J)*QA(J,NL1(J))/5760.
289    c      ENDDO      
290    C
291    c      IF(npas.eq.5760) then
292    c        open(73,file='transmoy',form='unformatted')
293    c        write(73) clsum
294    c        close(73)
295    c      ENDIF
296    Cchdbg
297    C      
298    C *********************************************************************
299    C *********************************************************************
300    
301          ABWLW1=0.7
302          DO 516 J=1,NGP
303           IF ( NLEVxy(J) .GT. 0 ) THEN
304            TAU(J,NLEVxy(J))=EXP(-(ABSLW+ABWLW*QA(J,NLEVxy(J)))*PSA(J)
305    cchdbg        TAU(J,NLEVxy(J))=EXP(-(ABSLW+ABWLW1*QA(J,NLEVxy(J)))*PSA(J)
306         *                     *DSIG(NLEVxy(J)))
307           ENDIF
308      516 CONTINUE
309    
310    C
311    C---
312          RETURN
313          END
314    
315    
316          SUBROUTINE RADLW (IMODE,TA,TS,ST4S,
317         *                  FTOP,FSFC,DFABS,FDOWN)
318    C--
319    C--   SUBROUTINE RADLW (IMODE,TA,TS,ST4S,
320    C--  *                  FTOP,FSFC,DFABS)
321    C--
322    C--   Purpose: Compute the absorption of longwave radiation
323    C--   Input:   IMODE  = index for operation mode (see below)
324    C--            TA     = absolute temperature (3-dim)
325    C--            TS     = surface temperature (2-dim)        [if IMODE=1]
326    C--            ST4S   = surface blackbody emission (2-dim) [if IMODE=2]
327    C--   Output:  ST4S   = surface blackbody emission (2-dim) [if IMODE=1]
328    C--            FTOP   = outgoing flux of lw rad. at the atm. top    (2-dim)
329    C--            FSFC   = net upw. flux of lw rad. at the surface     (2-dim)
330    C--            DFABS  = flux of lw rad. absorbed by each atm. layer (3-dim)
331    C--            FDOWN  = downward flux of lw rad. at the surface     (2-dim)
332    C--
333    
334    
335          IMPLICIT rEAL*8 (A-H,O-Z)
336    
337    
338    C     Resolution parameters
339    C
340    #include "atparam.h"
341    #include "atparam1.h"
342    #include "Lev_def.h"
343    C
344          PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT )
345    C
346    C     Constants + functions of sigma and latitude
347    C
348    #include "com_physcon.h"
349    C
350    C     Radiation parameters
351    C
352    #include "com_radcon.h"
353    C
354          REAL TA(NGP,NLEV), TS(NGP), ST4S(NGP)
355    
356          REAL FTOP(NGP), FSFC(NGP), DFABS(NGP,NLEV)
357          REAL FDOWN(NGP)
358    
359          REAL FLUX(NGP), BRAD(NGP), STCOR(NGP)
360          INTEGER NL1(NGP)
361    C
362    Cchdbg
363          INteger npas
364          SAVE npas
365          LOGICAL Ifirst
366          SAVE IFIRST
367          DATA Ifirst/.TRUE./
368          REAL FluxMoy(NGP)
369          REAL ST4SMoy(NGP)
370          SAVE FluxMoy, ST4SMoy
371    Cchdbg
372          DO J=1,NGP
373            NL1(J)=NLEVxy(J)-1
374          ENDDO
375    C
376          DO K=1,NLEV
377            DO J=1,NGP
378              DFABS(J,K)=0.
379            ENDDO
380          ENDDO
381    C
382    C---  1. Blackbody emission from atmospheric full and half levels.
383    C        Temperature is interpolated as a linear function of ln sigma.
384    C        At the lower boundary, the emission is linearly extrapolated;
385    C        at the upper boundary, the atmosphere is assumed isothermal.
386    C
387          DO 102 J=1,NGP
388          DO 102 K=1,NLEVxy(J)
389            ST4A(J,K,1)=TA(J,K)*TA(J,K)
390            ST4A(J,K,1)=SBC*ST4A(J,K,1)*ST4A(J,K,1)
391      102 CONTINUE
392    C
393          DO 104 J=1,NGP
394          DO 104 K=1,NL1(J)
395            ST4A(J,K,2)=TA(J,K)+WVI(K,2)*(TA(J,K+1)-TA(J,K))
396            ST4A(J,K,2)=ST4A(J,K,2)*ST4A(J,K,2)
397            ST4A(J,K,2)=SBC*ST4A(J,K,2)*ST4A(J,K,2)
398      104 CONTINUE
399    C
400          DO 106 J=1,NGP
401           IF ( NLEVxy(J) .GT. 0 ) THEN
402            ST4A(J,NLEVxy(J),2)=2.*ST4A(J,NLEVxy(J),1)-ST4A(J,NL1(J),2)
403           ENDIF
404      106 CONTINUE
405    C
406    C---  2. Empirical stratospheric correction
407    C
408          COR0= -13.
409          COR1=   0.
410          COR2=  24.
411    C
412          J0=0
413          DO JL=1,NLAT
414            CORR=COR0+COR1*FMU(JL,1)+COR2*FMU(JL,2)
415            DO J=J0+1,J0+NLON
416              STCOR(J)=CORR
417            ENDDO
418            J0=J0+NLON
419          ENDDO
420    C
421    C---  3. Emission ad absorption of longwave downward flux.
422    C        Downward emission is an average of the emission from the full level
423    C        and the half-level below, weighted according to the transmissivity
424    C        of the layer.
425    C      
426    C     3.1  Stratosphere
427    C
428          DO 312 J=1,NGP
429            BRAD(J)=ST4A(J,1,2)+TAU(J,1)*(ST4A(J,1,1)-ST4A(J,1,2))
430            FLUX(J)=(1.-TAU(J,1))*BRAD(J)
431            DFABS(J,1)=STCOR(J)-FLUX(J)
432      312 CONTINUE
433    C      
434    C     3.2  Troposphere
435    C
436          DO 322 J=1,NGP
437          DO 322 K=2,NLEVxy(J)
438            DFABS(J,K)=FLUX(J)
439            BRAD(J)=ST4A(J,K,2)+TAU(J,K)*(ST4A(J,K,1)-ST4A(J,K,2))
440            FLUX(J)=TAU(J,K)*(FLUX(J)-BRAD(J))+BRAD(J)
441            DFABS(J,K)=DFABS(J,K)-FLUX(J)
442      322 CONTINUE
443    C
444    C---  4. Emission ad absorption of longwave upward flux
445    C        Upward emission is an average of the emission from the full level
446    C        and the half-level above, weighted according to the transmissivity
447    C        of the layer (for the top layer, full-level emission is used).
448    C        Surface lw emission in the IR window goes directly into FTOP.
449    C      
450    C     4.1  Surface
451    C
452          IF (IMODE.LE.1) THEN
453            DO 412 J=1,NGP
454              ST4S(J)=TS(J)*TS(J)
455              ST4S(J)=SBC*ST4S(J)*ST4S(J)
456      412   CONTINUE
457          ENDIF
458    C
459    C **************************************************************
460    Cchdbg
461          if(ifirst) then
462            DO J=1,NGP
463              ST4SMoy(J)=0.
464              FluxMoy(J)=0.
465            ENDDO
466            npas=0.
467            ifirst=.FALSE.
468          endif
469    C
470          npas=npas+1
471          DO 413 J=1,NGP
472            ST4SMoy(J)=ST4SMoy(J)+ ST4S(J)
473            FluxMoy(J)=FluxMoy(J)+ Flux(J)
474      413 CONTINUE
475    C
476          if(npas.eq.5760) then
477            DO J=1,NGP
478              ST4SMoy(J)=ST4SMoy(J)/float(npas)
479              FluxMoy(J)=FluxMoy(J)/float(npas)
480            ENDDO
481            open(73,file='ST4Smoy',form='unformatted')
482            write(73) ST4SMoy
483            close(73)
484            open(74,file='FluxMoy',form='unformatted')
485            write(74) FluxMoy
486            close(74)
487          ENDIF
488    Cchdbg
489    C ****************************************************************
490    C
491    C
492          DO 414 J=1,NGP
493            FSFC(J)=ST4S(J)-FLUX(J)
494            FDOWN(J)=FLUX(J)
495            FTOP(J)=EPSLW*ST4S(J)
496            FLUX(J)=ST4S(J)-FTOP(J)
497      414 CONTINUE
498    C      
499    C     4.2  Troposphere
500    C
501          DO 422 J=1,NGP
502          DO 422 K=NLEVxy(J),2,-1
503            DFABS(J,K)=DFABS(J,K)+FLUX(J)
504            BRAD(J)=ST4A(J,K-1,2)+TAU(J,K)*(ST4A(J,K,1)-ST4A(J,K-1,2))
505            FLUX(J)=TAU(J,K)*(FLUX(J)-BRAD(J))+BRAD(J)
506            DFABS(J,K)=DFABS(J,K)-FLUX(J)
507      422 CONTINUE
508    C      
509    C     4.3  Stratosphere
510    C
511          DO 432 J=1,NGP
512            DFABS(J,1)=DFABS(J,1)+FLUX(J)
513            FLUX(J)=TAU(J,1)*(FLUX(J)-ST4A(J,1,1))+ST4A(J,1,1)
514            DFABS(J,1)=DFABS(J,1)-FLUX(J)
515      432 CONTINUE
516    C
517    C     4.4  Outgoing longwave radiation
518    C
519          DO 442 J=1,NGP
520    cdj     FTOP(J)=FTOP(J)+FLUX(J)
521            FTOP(J)=FTOP(J)+FLUX(J)-STCOR(J)
522      442 CONTINUE
523    cdj
524    c     write(0,*)'position j=20'
525    c     j=20
526    c     write(0,*)'ftop fsfc ftop-fsfc'
527    c     write(0,*)ftop(j),fsfc(j),ftop(j)-fsfc(j)
528    c     write(0,*)
529    c     write(0,*)'k dfabs'
530    c     do k = 1, nlevxy(j)
531    c       write(0,*)k,dfabs(j,k)
532    c     enddo
533    c     write(0,*)'sum dfabs'
534    c     write(0,*)sum(dfabs(j,:))
535    c     open(74,file='ftop0',form='unformatted',status='unknown')
536    c     write(74) ftop
537    c     open(75,file='stcor',form='unformatted',status='unknown')
538    c     write(75) stcor
539    c     stop
540    cdj
541    C
542    C---                                            
543          RETURN
544          END

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2

  ViewVC Help
Powered by ViewVC 1.1.22