source: trunk/LMDZ.MARS/libf/phymars/lwflux.F @ 3807

Last change on this file since 3807 was 3757, checked in by emillour, 6 weeks ago

Mars PCM:
More code tidying: puting routines in modules and modernizing some old
constructs. Tested to not change results with respect to previous version.
EM

File size: 11.3 KB
Line 
1      module lwflux_mod
2     
3      implicit none
4     
5      contains
6     
7      subroutine lwflux (ig0,kdlon,kflev,dp
8     .                   ,bsurf,btop,blev,blay,dbsublay
9     .                   ,tlay, tlev, dt0      ! pour sortie dans g2d uniquement
10     .                   ,emis
11     .                   , tautotal,omegtotal,gtotal
12     .                   ,coolrate,fluxground,fluxtop
13     .                   ,netrad)
14
15c----------------------------------------------------------------------
16c     LWFLUX     computes the fluxes
17c----------------------------------------------------------------------
18
19      use dimradmars_mod, only: ndlo2, nir, ndlon, nuco2, nflev
20      use yomlw_h, only: nlaylte, xi, xi_ground, gcp
21      use lwdiff_mod, only: lwdiff
22      implicit none
23 
24c----------------------------------------------------------------------
25c         0.1   arguments
26c               ---------
27c                                                            inputs:
28c                                                            -------
29      integer,intent(in) :: ig0
30      integer,intent(in) :: kdlon           ! part of ngrid
31      integer,intent(in) :: kflev           ! part of nlayer
32
33      real,intent(in) :: dp (ndlo2,kflev)   ! layer pressure thickness (Pa)
34
35      real,intent(in) :: bsurf (ndlo2,nir) ! surface spectral planck function
36      real,intent(in) :: blev (ndlo2,nir,kflev+1) ! level   spectral planck function
37      real,intent(in) :: blay (ndlo2,nir,kflev) ! layer   spectral planck function
38      real,intent(in) :: btop (ndlo2,nir) ! top spectral planck function
39      real,intent(in) :: dbsublay (ndlo2,nir,2*kflev)  ! layer gradient spectral planck
40                                                       ! function in sub layers
41
42      real,intent(in) :: dt0 (ndlo2) ! surface temperature discontinuity
43      real,intent(in) :: tlay (ndlo2,kflev) ! layer temperature
44      real,intent(in) :: tlev (ndlo2,kflev+1) ! level temperature
45
46      real,intent(in) :: emis (ndlo2) ! surface emissivity
47
48      real,intent(in) :: tautotal(ndlo2,kflev,nir)  ! \   Total single scattering
49      real,intent(in) :: omegtotal(ndlo2,kflev,nir) !  >  properties (Addition of the
50      real,intent(in) :: gtotal(ndlo2,kflev,nir)    ! /   NAERKIND aerosols prop.)
51
52
53c                                                            outputs:
54c                                                            --------
55      real,intent(out) :: coolrate(ndlo2,kflev) ! radiative cooling rate (K/s)
56      real,intent(out) :: netrad (ndlo2,kflev) ! radiative budget (W/m2)
57      real,intent(out) :: fluxground(ndlo2) ! downward flux on the ground
58                                            ! for surface radiative budget
59      real,intent(out) :: fluxtop(ndlo2) ! upward flux on the top of atm ("OLR")
60
61
62c----------------------------------------------------------------------
63c         0.2   local arrays
64c               ------------
65
66      integer ja,jl,j,i,ig1d,ig,l
67      real  ksidb (ndlon,nuco2+1,0:nflev+1,0:nflev+1) ! net exchange rate (W/m2)
68
69      real dpsgcp (0:nflev+1,0:nflev+1)    !  dp/(g.cp)
70      real temp (0:nflev+1,0:nflev+1)       
71
72      real fluxdiff(ndlon,2,nflev+1)  ! diffusion flux: upward(1) downward(2)
73
74      real*4 reel4
75
76c     To compute IR flux in the atmosphere  (For diagnostic only !!)
77      logical computeflux
78      real coefd(kdlon,nuco2,nflev+1,nflev+1)
79      real coefu(kdlon,nuco2,0:nflev,nflev+1)
80      real flw_up(kdlon,nflev+1), flw_dn(kdlon,nflev+1) ! fluxes (W/m2)
81
82
83      ksidb(:,:,:,:)=0
84
85c----------------------------------------------------------------------
86c         1.1   exchanges (layer i <--> all layers up to i)
87c               -------------------------------------------
88
89      do i = 1,nlaylte
90        do j = i+1,nlaylte
91          do ja = 1,nuco2
92            do jl = 1,kdlon
93
94      ksidb(jl,ja,i,j) = xi(ig0+jl,ja,i,j)
95     .                 * (blay(jl,ja,j)-blay(jl,ja,i))
96c                                                        ksidb reciprocity
97c                                                        -----------------
98      ksidb(jl,ja,j,i) = -ksidb(jl,ja,i,j)
99     
100            enddo
101          enddo
102        enddo
103      enddo
104
105c----------------------------------------------------------------------
106c         1.2   exchanges (ground <--> all layers)
107c               ----------------------------------
108
109      do i = 1,nlaylte
110        do ja = 1,nuco2
111          do jl = 1,kdlon
112
113      ksidb(jl,ja,i,0) = xi(ig0+jl,ja,0,i)
114     .                 * (bsurf(jl,ja)-blay(jl,ja,i))
115c                                                        ksidb reciprocity
116c                                                        -----------------
117      ksidb(jl,ja,0,i) = -ksidb(jl,ja,i,0)
118
119          enddo
120        enddo
121      enddo
122
123c--------------------------------------------------------
124c  Here we add the neighbour contribution
125c           for exchanges between ground and first layer
126c--------------------------------------------------------
127
128        do ja = 1,nuco2
129          do jl = 1,kdlon
130
131      ksidb(jl,ja,1,0) = ksidb(jl,ja,1,0)
132     .                 - xi_ground(ig0+jl,ja)
133     .                 * (blev(jl,ja,1)-blay(jl,ja,1))
134
135cc                                                       ksidb reciprocity
136cc                                                       -----------------
137      ksidb(jl,ja,0,1) = - ksidb(jl,ja,1,0)
138
139          enddo
140        enddo
141
142c----------------------------------------------------------------------
143c         1.3   exchanges (layer i <--> space)
144c               ------------------------------
145
146      do i = 1,nlaylte
147        do ja = 1,nuco2
148          do jl = 1,kdlon
149
150      ksidb(jl,ja,i,nlaylte+1) = xi(ig0+jl,ja,i,nlaylte+1)
151     .                       * (-blay(jl,ja,i))
152c                                                        ksidb reciprocity
153c                                                        -----------------
154      ksidb(jl,ja,nlaylte+1,i) = - ksidb(jl,ja,i,nlaylte+1)
155
156          enddo
157        enddo
158      enddo
159
160c----------------------------------------------------------------------
161c         1.4   exchanges (ground <--> space)
162c               -----------------------------
163
164      do ja = 1,nuco2
165        do jl = 1,kdlon
166
167      ksidb(jl,ja,0,nlaylte+1) = xi(ig0+jl,ja,0,nlaylte+1)
168     .                       * (-bsurf(jl,ja))
169
170c                                                        ksidb reciprocity
171c                                                        -----------------
172      ksidb(jl,ja,nlaylte+1,0) = - ksidb(jl,ja,0,nlaylte+1)
173
174        enddo
175      enddo
176
177c----------------------------------------------------------------------
178c         2.0   sum of band 1 and 2 of co2 contribution
179c               ---------------------------------------
180
181      do i = 0,nlaylte+1
182        do j = 0,nlaylte+1
183          do jl = 1,kdlon
184
185            ksidb(jl,3,i,j)= ksidb(jl,1,i,j) + ksidb(jl,2,i,j)
186
187          enddo
188        enddo
189      enddo
190
191c----------------------------------------------------------------------
192c         3.0   Diffusion
193c               ---------
194
195      i = nlaylte+1
196      do jl = 1,kdlon
197        fluxdiff(jl,1,i)   = 0.
198        fluxdiff(jl,2,i)   = 0.
199      enddo
200
201      call lwdiff (kdlon,kflev
202     .         ,bsurf,btop,dbsublay
203     .         ,tautotal,omegtotal,gtotal
204     .         ,emis,fluxdiff)
205
206c----------------------------------------------------------------------
207c         4.0   Radiative Budget for each layer i
208c               ---------------------------------
209
210      do i = 1,nlaylte
211        do jl = 1,kdlon
212            netrad(jl,i) = 0.
213        enddo
214      enddo
215
216      do i = 1,nlaylte
217        do j = 0,nlaylte+1
218          do jl = 1,kdlon
219
220            netrad(jl,i) = netrad(jl,i) + ksidb(jl,3,i,j)
221
222          enddo
223        enddo
224      enddo
225c                                                 diffusion contribution
226c                                                 ----------------------
227      do i = 1,nlaylte
228        do jl = 1,kdlon
229
230          netrad(jl,i) = netrad(jl,i)
231     .                 - fluxdiff(jl,1,i+1) - fluxdiff(jl,2,i+1)   
232     .                 + fluxdiff(jl,1,i)   + fluxdiff(jl,2,i)
233
234        enddo
235      enddo
236
237c----------------------------------------------------------------------
238c         4.0   cooling rate for each layer i
239c               -----------------------------
240
241      do i = 1,nlaylte
242        do jl = 1,kdlon
243
244          coolrate(jl,i) = gcp * netrad(jl,i) / dp(jl,i)
245
246        enddo
247      enddo
248
249c----------------------------------------------------------------------
250c         5.0   downward flux (all layers --> ground): "fluxground"
251c               ---------------------------------------------------
252
253      do jl = 1,kdlon
254          fluxground(jl) = 0.
255      enddo
256
257      do i = 1,nlaylte
258        do ja = 1,nuco2
259          do jl = 1,kdlon
260
261      fluxground(jl) = fluxground(jl)
262     .               + xi(ig0+jl,ja,0,i) * (blay(jl,ja,i))
263
264          enddo
265        enddo
266      enddo
267   
268      do jl = 1,kdlon
269        fluxground(jl) = fluxground(jl) - fluxdiff(jl,2,1) 
270      enddo
271
272c----------------------------------------------------------------------
273c         6.0   outgoing flux (all layers --> space): "fluxtop"
274c               ---------------------------------------------------
275
276      do jl = 1,kdlon
277          fluxtop(jl) = 0.
278      enddo
279
280      do i = 0,nlaylte
281        do jl = 1,kdlon
282           fluxtop(jl) = fluxtop(jl)- ksidb(jl,3,i,nlaylte+1)
283        enddo
284      enddo
285   
286      do jl = 1,kdlon
287        fluxtop(jl) = fluxtop(jl) + fluxdiff(jl,1,nlaylte+1) 
288      enddo
289
290c----------------------------------------------------------------------
291c         6.5 ONLY FOR DIAGNOSTIC : Compute IR flux in the atmosphere
292c             -------------------
293c     The broadband fluxes (W.m-2) at every level  from surface level (l=1)
294c     up the top of the  upper layer (here: l=nlaylte+1) are:
295c     upward : flw_up(ig1d,l)   ;   downward : flw_dn(ig1d,j)
296c
297      computeflux = .false.
298
299      IF (computeflux) THEN  ! not used by the GCM only for diagnostic !
300c      upward flux
301c      ~~~~~~~~~~~
302       do i = 0,nlaylte
303        do j = 1,nlaylte+1
304         do ja = 1,nuco2
305          do jl = 1,kdlon
306            coefu(jl,ja,i,j) =0.
307            do l=j,nlaylte+1
308              coefu(jl,ja,i,j)=coefu(jl,ja,i,j)+xi(ig0+jl,ja,l,i)
309            end do
310
311          enddo
312         enddo
313        enddo
314       enddo
315       do j = 1,nlaylte+1
316        do jl = 1,kdlon
317           flw_up(jl,j) = 0.
318           do ja = 1,nuco2
319             flw_up(jl,j)=flw_up(jl,j)+bsurf(jl,ja)*coefu(jl,ja,0,j)
320             do i=1,j-1
321              flw_up(jl,j)=flw_up(jl,j)+blay(jl,ja,i)*coefu(jl,ja,i,j)
322             end do
323           end do
324           flw_up(jl,j)=flw_up(jl,j) + fluxdiff(jl,1,j)
325        end do
326       end do
327
328c      downward flux
329c      ~~~~~~~~~~~~~
330       do i = 1,nlaylte+1
331        do j = 1,nlaylte+1
332         do ja = 1,nuco2
333          do jl = 1,kdlon
334            coefd(jl,ja,i,j) =0.
335            do l=0,j-1
336              coefd(jl,ja,i,j)=coefd(jl,ja,i,j)+xi(ig0+jl,ja,l,i)
337            end do
338          enddo
339         enddo
340        enddo
341       enddo
342       do j = 1,nlaylte+1
343        do jl = 1,kdlon
344           flw_dn(jl,j) = 0.
345           do ja = 1,nuco2
346             do i=j,nlaylte
347              flw_dn(jl,j)=flw_dn(jl,j)+blay(jl,ja,i)*coefd(jl,ja,i,j)
348             end do
349           end do
350           flw_dn(jl,j)=flw_dn(jl,j) - fluxdiff(jl,2,j)
351        end do
352       end do
353      END IF ! of IF (computeflux)
354
355      end subroutine lwflux
356     
357      end module lwflux_mod
Note: See TracBrowser for help on using the repository browser.