source: trunk/LMDZ.GENERIC/libf/phystd/gfluxi.F @ 997

Last change on this file since 997 was 995, checked in by aslmd, 11 years ago

LMDZ.GENERIC. Updated saturn 1D files. Added comments around the gluxi problem (will be further investigated). removed abort because Helium rayleigh scattering is missing

File size: 7.5 KB
RevLine 
[135]1      SUBROUTINE GFLUXI(NLL,TLEV,NW,DW,DTAU,TAUCUM,W0,COSBAR,UBARI,
2     *                  RSF,BTOP,BSURF,FTOPUP,FMIDP,FMIDM)
3
4C  THIS SUBROUTINE TAKES THE OPTICAL CONSTANTS AND BOUNDARY CONDITIONS
5C  FOR THE INFRARED FLUX AT ONE WAVELENGTH AND SOLVES FOR THE FLUXES AT
6C  THE LEVELS. THIS VERSION IS SET UP TO WORK WITH LAYER OPTICAL DEPTHS
7C  MEASURED FROM THE TOP OF EACH LAYER.  THE TOP OF EACH LAYER HAS 
8C  OPTICAL DEPTH ZERO.  IN THIS SUB LEVEL N IS ABOVE LAYER N. THAT IS LAYER N
9C  HAS LEVEL N ON TOP AND LEVEL N+1 ON BOTTOM. OPTICAL DEPTH INCREASES
10C  FROM TOP TO BOTTOM. SEE C.P. MCKAY, TGM NOTES.
11C THE TRI-DIAGONAL MATRIX SOLVER IS DSOLVER AND IS DOUBLE PRECISION SO MANY
12C VARIABLES ARE PASSED AS SINGLE THEN BECOME DOUBLE IN DSOLVER
13C
14C NLL            = NUMBER OF LEVELS (NLAYERS + 1) MUST BE LESS THAT NL (101)
15C TLEV(L_LEVELS) = ARRAY OF TEMPERATURES AT GCM LEVELS
16C WAVEN          = WAVELENGTH FOR THE COMPUTATION
17C DW             = WAVENUMBER INTERVAL
18C DTAU(NLAYER)   = ARRAY OPTICAL DEPTH OF THE LAYERS
19C W0(NLEVEL)     = SINGLE SCATTERING ALBEDO
20C COSBAR(NLEVEL) = ASYMMETRY FACTORS, 0=ISOTROPIC
21C UBARI          = AVERAGE ANGLE, MUST BE EQUAL TO 0.5 IN IR
22C RSF            = SURFACE REFLECTANCE
23C BTOP           = UPPER BOUNDARY CONDITION ON IR INTENSITY (NOT FLUX)
24C BSURF          = SURFACE EMISSION = (1-RSFI)*PLANCK, INTENSITY (NOT FLUX)
25C FP(NLEVEL)     = UPWARD FLUX AT LEVELS
26C FM(NLEVEL)     = DOWNWARD FLUX AT LEVELS
27C FMIDP(NLAYER)  = UPWARD FLUX AT LAYER MIDPOINTS
28C FMIDM(NLAYER)  = DOWNWARD FLUX AT LAYER MIDPOINTS
29C
30C----------------------------------------------------------------------C
31
32      use radinc_h
33      use radcommon_h, only: planckir
34
35      implicit none
36
37#include "comcstfi.h"
38
39      INTEGER NLP
40      PARAMETER (NLP=101) ! MUST BE LARGER THAN NLEVEL
41
42      INTEGER NLL, NLAYER, L, NW, NT, NT2
43      REAL*8  TERM, CPMID, CMMID
44      REAL*8  PLANCK
45      REAL*8  EM,EP
46      REAL*8  COSBAR(L_NLAYRAD), W0(L_NLAYRAD), DTAU(L_NLAYRAD)
47      REAL*8  TAUCUM(L_LEVELS), DTAUK
48      REAL*8  TLEV(L_LEVELS)
49      REAL*8  WAVEN, DW, UBARI, RSF
50      REAL*8  BTOP, BSURF, FMIDP(L_NLAYRAD), FMIDM(L_NLAYRAD)
51      REAL*8  B0(NLP),B1(NLP),ALPHA(NLP),LAMDA(NLP),XK1(NLP),XK2(NLP)
52      REAL*8  GAMA(NLP),CP(NLP),CM(NLP),CPM1(NLP),CMM1(NLP),E1(NLP)
53      REAL*8  E2(NLP),E3(NLP),E4(NLP)
54
55      REAL*8  FTOPUP, FLUXUP, FLUXDN
56
57      real*8 :: TAUMAX = L_TAUMAX
58
59C======================================================================C
60 
61C     WE GO WITH THE HEMISPHERIC CONSTANT APPROACH IN THE INFRARED
62 
63
[543]64      IF (NLL .GT. NLP) STOP 'PARAMETER NL TOO SMALL IN GFLUXI'
[135]65
66      NLAYER = L_NLAYRAD
67
68      DO L=1,L_NLAYRAD-1
[804]69
70!----------------------------------------------------
71! There is a problem when W0 = 1
72!         open(888,file='W0')
73!           if ((W0(L).eq.0.).or.(W0(L).eq.1.)) then
74!             write(888,*) W0(L), L, 'gfluxi'
75!           endif
76! Prevent this with an if statement:
[959]77        if (W0(L).eq.1.D0) then
78           W0(L) = 0.99999D0
[804]79        endif
80!----------------------------------------------------
81
[959]82        ALPHA(L) = SQRT( (1.0D0-W0(L))/(1.0D0-W0(L)*COSBAR(L)) )
83        LAMDA(L) = ALPHA(L)*(1.0D0-W0(L)*COSBAR(L))/UBARI
[135]84
[543]85        NT    = int(TLEV(2*L)*NTfac)   - NTstar+1
86        NT2   = int(TLEV(2*L+2)*NTfac) - NTstar+1
[135]87
88        B1(L) = (PLANCKIR(NW,NT2)-PLANCKIR(NW,NT))/DTAU(L)
89        B0(L) = PLANCKIR(NW,NT)
90      END DO
91
92C     Take care of special lower layer
93
94      L        = L_NLAYRAD
[804]95
96      if (W0(L).eq.1.) then
[959]97          W0(L) = 0.99999D0
[804]98      end if
99
[959]100      ALPHA(L) = SQRT( (1.0D0-W0(L))/(1.0D0-W0(L)*COSBAR(L)) )
101      LAMDA(L) = ALPHA(L)*(1.0D0-W0(L)*COSBAR(L))/UBARI
[135]102
[995]103      ! Tsurf is used for 1st layer source function
104      ! -- same results for most thin atmospheres
105      ! -- and stabilizes integrations
[543]106      NT    = int(TLEV(2*L+1)*NTfac) - NTstar+1
[995]107
108      !! For deep, opaque, thick first layers (e.g. Saturn)
109      !! what is below works much better, not unstable, ...
110      !! ... and actually fully accurate because 1st layer temp (JL)
111      !NT    = int(TLEV(2*L)*NTfac) - NTstar+1
112      !! (or this one yields same results
113      !NT    = int( (TLEV(2*L)+TLEV(2*L+1))*0.5*NTfac ) - NTstar+1
114
[543]115      NT2   = int(TLEV(2*L)*NTfac)   - NTstar+1
[135]116      B1(L) = (PLANCKIR(NW,NT)-PLANCKIR(NW,NT2))/DTAU(L)
117      B0(L) = PLANCKIR(NW,NT2)
118
119      DO L=1,L_NLAYRAD
[959]120        GAMA(L) = (1.0D0-ALPHA(L))/(1.0D0+ALPHA(L))
121        TERM    = UBARI/(1.0D0-W0(L)*COSBAR(L))
[135]122
123C       CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
124C       BOTTOM OF THE LAYER.  THAT IS AT DTAU OPTICAL DEPTH
125
126        CP(L) = B0(L)+B1(L)*DTAU(L) +B1(L)*TERM
127        CM(L) = B0(L)+B1(L)*DTAU(L) -B1(L)*TERM
128
129C       CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED
130C       AT THE TOP OF THE LAYER, THAT IS ZERO OPTICAL DEPTH
131
132        CPM1(L) = B0(L)+B1(L)*TERM
133        CMM1(L) = B0(L)-B1(L)*TERM
134      END DO
135 
136C     NOW CALCULATE THE EXPONENTIAL TERMS NEEDED
137C     FOR THE TRIDIAGONAL ROTATED LAYERED METHOD
138C     WARNING IF DTAU(J) IS GREATER THAN ABOUT 35 (VAX)
139C     WE CLIP IT TO AVOID OVERFLOW.
140
141      DO L=1,L_NLAYRAD
142
143C       CLIP THE EXPONENTIAL HERE.
144
145        EP    = EXP( MIN((LAMDA(L)*DTAU(L)),TAUMAX))
[959]146        EM    = 1.0D0/EP
[135]147        E1(L) = EP+GAMA(L)*EM
148        E2(L) = EP-GAMA(L)*EM
149        E3(L) = GAMA(L)*EP+EM
150        E4(L) = GAMA(L)*EP-EM
151      END DO
152 
153c     B81=BTOP  ! RENAME BEFORE CALLING DSOLVER - used to be to set
154c     B82=BSURF ! them to real*8 - but now everything is real*8
155c     R81=RSF   ! so this may not be necessary
156
157C     Double precision tridiagonal solver
158
159      CALL DSOLVER(NLAYER,GAMA,CP,CM,CPM1,CMM1,E1,E2,E3,E4,BTOP,
160     *             BSURF,RSF,XK1,XK2)
161 
162
163C     NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS.
164 
165      DO L=1,L_NLAYRAD-1
166        DTAUK = TAUCUM(2*L+1)-TAUCUM(2*L)
167        EP    = EXP(MIN(LAMDA(L)*DTAUK,TAUMAX)) ! CLIPPED EXPONENTIAL
[959]168        EM    = 1.0D0/EP
[961]169        TERM  = UBARI/(1.D0-W0(L)*COSBAR(L))
[135]170
171C       CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
172C       BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH
173
174        CPMID    = B0(L)+B1(L)*DTAUK +B1(L)*TERM
175        CMMID    = B0(L)+B1(L)*DTAUK -B1(L)*TERM
176        FMIDP(L) = XK1(L)*EP + GAMA(L)*XK2(L)*EM + CPMID
177        FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID
178
179C       FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
180
181        FMIDP(L) = FMIDP(L)*PI
182        FMIDM(L) = FMIDM(L)*PI
183      END DO
184 
185C     And now, for the special bottom layer
186
187      L    = L_NLAYRAD
188
189      EP   = EXP(MIN((LAMDA(L)*DTAU(L)),TAUMAX)) ! CLIPPED EXPONENTIAL
[959]190      EM   = 1.0D0/EP
191      TERM = UBARI/(1.D0-W0(L)*COSBAR(L))
[135]192
193C     CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
194C     BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH
195
196      CPMID    = B0(L)+B1(L)*DTAU(L) +B1(L)*TERM
197      CMMID    = B0(L)+B1(L)*DTAU(L) -B1(L)*TERM
198      FMIDP(L) = XK1(L)*EP + GAMA(L)*XK2(L)*EM + CPMID
199      FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID
200 
201C     FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
202
203      FMIDP(L) = FMIDP(L)*PI
204      FMIDM(L) = FMIDM(L)*PI
205
206C     FLUX AT THE PTOP LEVEL
207
[959]208      EP   = 1.0D0
209      EM   = 1.0D0
210      TERM = UBARI/(1.0D0-W0(1)*COSBAR(1))
[135]211
212C     CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
213C     BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH
214
215      CPMID  = B0(1)+B1(1)*TERM
216      CMMID  = B0(1)-B1(1)*TERM
217
218      FLUXUP = XK1(1)*EP + GAMA(1)*XK2(1)*EM + CPMID
219      FLUXDN = XK1(1)*EP*GAMA(1) + XK2(1)*EM + CMMID
220 
221C     FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT
222
223      FTOPUP = (FLUXUP-FLUXDN)*PI
224
225
226      RETURN
227      END
Note: See TracBrowser for help on using the repository browser.