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

Last change on this file since 1242 was 1047, checked in by emillour, 11 years ago

Mars GCM:

  • IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
  • Also updated EOF routines, everything is now in eofdump_mod.F90
  • Removed unused routine lectfux.F (in dyn3d)

EM

File size: 15.4 KB
Line 
1       subroutine lwflux (ig0,kdlon,kflev,dp
2     .                   ,bsurf,btop,blev,blay,dbsublay
3     .                   ,tlay, tlev, dt0      ! pour sortie dans g2d uniquement
4     .                   ,emis
5     .                   , tautotal,omegtotal,gtotal
6     .                   ,coolrate,fluxground,fluxtop
7     .                   ,netrad)
8
9c----------------------------------------------------------------------
10c     LWFLUX     computes the fluxes
11c----------------------------------------------------------------------
12
13      use dimradmars_mod, only: ndlo2, nir, ndlon, nuco2, nflev
14      use yomlw_h, only: nlaylte, xi, xi_ground, gcp
15      implicit none
16 
17
18!#include "dimensions.h"
19!#include "dimphys.h"
20!#include "dimradmars.h"
21#include "callkeys.h"
22#include "comg1d.h"
23 
24!#include "yomlw.h"
25
26c----------------------------------------------------------------------
27c         0.1   arguments
28c               ---------
29c                                                            inputs:
30c                                                            -------
31      integer ig0
32      integer kdlon                 ! part of ngrid
33      integer kflev                 ! part of nlayer
34
35      real dp (ndlo2,kflev)         ! layer pressure thickness (Pa)
36
37      real bsurf (ndlo2,nir)             ! surface spectral planck function
38      real blev (ndlo2,nir,kflev+1)      ! level   spectral planck function
39      real blay (ndlo2,nir,kflev)        ! layer   spectral planck function
40      real btop (ndlo2,nir)              ! top spectral planck function
41      real dbsublay (ndlo2,nir,2*kflev)  ! layer gradient spectral planck
42                                         ! function in sub layers
43
44      real dt0 (ndlo2)                 ! surface temperature discontinuity
45      real tlay (ndlo2,kflev)          ! layer temperature
46      real tlev (ndlo2,kflev+1)        ! level temperature
47
48      real emis (ndlo2)                  ! surface emissivity
49
50      real tautotal(ndlo2,kflev,nir)  ! \   Total single scattering
51      real omegtotal(ndlo2,kflev,nir) !  >  properties (Addition of the
52      real gtotal(ndlo2,kflev,nir)    ! /   NAERKIND aerosols prop.)
53
54
55c                                                            outputs:
56c                                                            --------
57      real coolrate(ndlo2,kflev)    ! radiative cooling rate (K/s)
58      real netrad (ndlo2,kflev)     ! radiative budget (W/m2)
59      real fluxground(ndlo2)        ! downward flux on the ground
60                                    ! for surface radiative budget
61      real fluxtop(ndlo2)           ! upward flux on the top of atm ("OLR")
62
63
64c----------------------------------------------------------------------
65c         0.2   local arrays
66c               ------------
67
68      integer ja,jl,j,i,ig1d,ig,l,ndim
69!      parameter(ndim = ndlon*(nuco2+1)*(nflev+2)*(nflev+2))
70      real  ksidb (ndlon,nuco2+1,0:nflev+1,0:nflev+1) ! net exchange rate (W/m2)
71
72      real dpsgcp (0:nflev+1,0:nflev+1)    !  dp/(g.cp)
73      real temp (0:nflev+1,0:nflev+1)       
74
75      real fluxdiff(ndlon,2,nflev+1)  ! diffusion flux: upward(1) downward(2)
76
77      real*4 reel4
78
79c     To compute IR flux in the atmosphere  (For diagnostic only !!)
80      logical computeflux
81      real coefd(kdlon,nuco2,nflev+1,nflev+1)
82      real coefu(kdlon,nuco2,0:nflev,nflev+1)
83      real flw_up(kdlon,nflev+1), flw_dn(kdlon,nflev+1) ! fluxes (W/m2)
84
85
86      ndim = ndlon*(nuco2+1)*(nflev+2)*(nflev+2)
87      call zerophys(ndim, ksidb)
88
89c----------------------------------------------------------------------
90c         1.1   exchanges (layer i <--> all layers up to i)
91c               -------------------------------------------
92
93      do i = 1,nlaylte
94        do j = i+1,nlaylte
95          do ja = 1,nuco2
96            do jl = 1,kdlon
97
98      ksidb(jl,ja,i,j) = xi(ig0+jl,ja,i,j)
99     .                 * (blay(jl,ja,j)-blay(jl,ja,i))
100c                                                        ksidb reciprocity
101c                                                        -----------------
102      ksidb(jl,ja,j,i) = -ksidb(jl,ja,i,j)
103     
104            enddo
105          enddo
106        enddo
107      enddo
108
109c----------------------------------------------------------------------
110c         1.2   exchanges (ground <--> all layers)
111c               ----------------------------------
112
113      do i = 1,nlaylte
114        do ja = 1,nuco2
115          do jl = 1,kdlon
116
117      ksidb(jl,ja,i,0) = xi(ig0+jl,ja,0,i)
118     .                 * (bsurf(jl,ja)-blay(jl,ja,i))
119c                                                        ksidb reciprocity
120c                                                        -----------------
121      ksidb(jl,ja,0,i) = -ksidb(jl,ja,i,0)
122
123          enddo
124        enddo
125      enddo
126
127c--------------------------------------------------------
128c  Here we add the neighbour contribution
129c           for exchanges between ground and first layer
130c--------------------------------------------------------
131
132        do ja = 1,nuco2
133          do jl = 1,kdlon
134
135      ksidb(jl,ja,1,0) = ksidb(jl,ja,1,0)
136     .                 - xi_ground(ig0+jl,ja)
137     .                 * (blev(jl,ja,1)-blay(jl,ja,1))
138
139cc                                                       ksidb reciprocity
140cc                                                       -----------------
141      ksidb(jl,ja,0,1) = - ksidb(jl,ja,1,0)
142
143          enddo
144        enddo
145
146c----------------------------------------------------------------------
147c         1.3   exchanges (layer i <--> space)
148c               ------------------------------
149
150      do i = 1,nlaylte
151        do ja = 1,nuco2
152          do jl = 1,kdlon
153
154      ksidb(jl,ja,i,nlaylte+1) = xi(ig0+jl,ja,i,nlaylte+1)
155     .                       * (-blay(jl,ja,i))
156c                                                        ksidb reciprocity
157c                                                        -----------------
158      ksidb(jl,ja,nlaylte+1,i) = - ksidb(jl,ja,i,nlaylte+1)
159
160          enddo
161        enddo
162      enddo
163
164c----------------------------------------------------------------------
165c         1.4   exchanges (ground <--> space)
166c               -----------------------------
167
168      do ja = 1,nuco2
169        do jl = 1,kdlon
170
171      ksidb(jl,ja,0,nlaylte+1) = xi(ig0+jl,ja,0,nlaylte+1)
172     .                       * (-bsurf(jl,ja))
173
174c                                                        ksidb reciprocity
175c                                                        -----------------
176      ksidb(jl,ja,nlaylte+1,0) = - ksidb(jl,ja,0,nlaylte+1)
177
178        enddo
179      enddo
180
181c----------------------------------------------------------------------
182c         2.0   sum of band 1 and 2 of co2 contribution
183c               ---------------------------------------
184
185      do i = 0,nlaylte+1
186        do j = 0,nlaylte+1
187          do jl = 1,kdlon
188
189            ksidb(jl,3,i,j)= ksidb(jl,1,i,j) + ksidb(jl,2,i,j)
190
191          enddo
192        enddo
193      enddo
194
195c----------------------------------------------------------------------
196c         3.0   Diffusion
197c               ---------
198
199      i = nlaylte+1
200      do jl = 1,kdlon
201        fluxdiff(jl,1,i)   = 0.
202        fluxdiff(jl,2,i)   = 0.
203      enddo
204
205      call lwdiff (kdlon,kflev
206     .         ,bsurf,btop,dbsublay
207     .         ,tautotal,omegtotal,gtotal
208     .         ,emis,fluxdiff)
209
210c----------------------------------------------------------------------
211c         4.0   Radiative Budget for each layer i
212c               ---------------------------------
213
214      do i = 1,nlaylte
215        do jl = 1,kdlon
216            netrad(jl,i) = 0.
217        enddo
218      enddo
219
220      do i = 1,nlaylte
221        do j = 0,nlaylte+1
222          do jl = 1,kdlon
223
224            netrad(jl,i) = netrad(jl,i) + ksidb(jl,3,i,j)
225
226          enddo
227        enddo
228      enddo
229c                                                 diffusion contribution
230c                                                 ----------------------
231      do i = 1,nlaylte
232        do jl = 1,kdlon
233
234          netrad(jl,i) = netrad(jl,i)
235     .                 - fluxdiff(jl,1,i+1) - fluxdiff(jl,2,i+1)   
236     .                 + fluxdiff(jl,1,i)   + fluxdiff(jl,2,i)
237
238        enddo
239      enddo
240
241c----------------------------------------------------------------------
242c         4.0   cooling rate for each layer i
243c               -----------------------------
244
245      do i = 1,nlaylte
246        do jl = 1,kdlon
247
248          coolrate(jl,i) = gcp * netrad(jl,i) / dp(jl,i)
249
250        enddo
251      enddo
252
253c----------------------------------------------------------------------
254c         5.0   downward flux (all layers --> ground): "fluxground"
255c               ---------------------------------------------------
256
257      do jl = 1,kdlon
258          fluxground(jl) = 0.
259      enddo
260
261      do i = 1,nlaylte
262        do ja = 1,nuco2
263          do jl = 1,kdlon
264
265      fluxground(jl) = fluxground(jl)
266     .               + xi(ig0+jl,ja,0,i) * (blay(jl,ja,i))
267
268          enddo
269        enddo
270      enddo
271   
272      do jl = 1,kdlon
273        fluxground(jl) = fluxground(jl) - fluxdiff(jl,2,1) 
274      enddo
275
276c----------------------------------------------------------------------
277c         6.0   outgoing flux (all layers --> space): "fluxtop"
278c               ---------------------------------------------------
279
280      do jl = 1,kdlon
281          fluxtop(jl) = 0.
282      enddo
283
284      do i = 0,nlaylte
285        do jl = 1,kdlon
286           fluxtop(jl) = fluxtop(jl)- ksidb(jl,3,i,nlaylte+1)
287        enddo
288      enddo
289   
290      do jl = 1,kdlon
291        fluxtop(jl) = fluxtop(jl) + fluxdiff(jl,1,nlaylte+1) 
292      enddo
293
294c----------------------------------------------------------------------
295c         6.5 ONLY FOR DIAGNOSTIC : Compute IR flux in the atmosphere
296c             -------------------
297c     The broadband fluxes (W.m-2) at every level  from surface level (l=1)
298c     up the top of the  upper layer (here: l=nlaylte+1) are:
299c     upward : flw_up(ig1d,l)   ;   downward : flw_dn(ig1d,j)
300c
301      computeflux = .false.
302
303      IF (computeflux) THEN  ! not used by the GCM only for diagnostic !
304c      upward flux
305c      ~~~~~~~~~~~
306       do i = 0,nlaylte
307        do j = 1,nlaylte+1
308         do ja = 1,nuco2
309          do jl = 1,kdlon
310            coefu(jl,ja,i,j) =0.
311            do l=j,nlaylte+1
312              coefu(jl,ja,i,j)=coefu(jl,ja,i,j)+xi(ig0+jl,ja,l,i)
313            end do
314
315          enddo
316         enddo
317        enddo
318       enddo
319       do j = 1,nlaylte+1
320        do jl = 1,kdlon
321           flw_up(jl,j) = 0.
322           do ja = 1,nuco2
323             flw_up(jl,j)=flw_up(jl,j)+bsurf(jl,ja)*coefu(jl,ja,0,j)
324             do i=1,j-1
325              flw_up(jl,j)=flw_up(jl,j)+blay(jl,ja,i)*coefu(jl,ja,i,j)
326             end do
327           end do
328           flw_up(jl,j)=flw_up(jl,j) + fluxdiff(jl,1,j)
329        end do
330       end do
331
332c      downward flux
333c      ~~~~~~~~~~~~~
334       do i = 1,nlaylte+1
335        do j = 1,nlaylte+1
336         do ja = 1,nuco2
337          do jl = 1,kdlon
338            coefd(jl,ja,i,j) =0.
339            do l=0,j-1
340              coefd(jl,ja,i,j)=coefd(jl,ja,i,j)+xi(ig0+jl,ja,l,i)
341            end do
342          enddo
343         enddo
344        enddo
345       enddo
346       do j = 1,nlaylte+1
347        do jl = 1,kdlon
348           flw_dn(jl,j) = 0.
349           do ja = 1,nuco2
350             do i=j,nlaylte
351              flw_dn(jl,j)=flw_dn(jl,j)+blay(jl,ja,i)*coefd(jl,ja,i,j)
352             end do
353           end do
354           flw_dn(jl,j)=flw_dn(jl,j) - fluxdiff(jl,2,j)
355        end do
356       end do
357      END IF
358
359c----------------------------------------------------------------------
360c         7.0   outputs Grads 2D
361c               ----------------
362
363c ig1d: point de la grille physique ou on veut faire la sortie
364c ig0+1:  point du decoupage de la grille physique
365
366c#ifdef undim
367      if (callg2d) then
368
369      ig1d = kdlon/2 + 1
370c     ig1d = kdlon
371
372      if ((ig0+1).LE.ig1d .and. ig1d.LE.(ig0+kdlon)
373     .    .OR.  kdlon.EQ.1   ) then
374
375          ig = ig1d-ig0
376        print*, 'Sortie g2d: ig1d, ig, ig0', ig1d, ig, ig0
377
378c--------------------------------------------
379c   Ouverture de g2d.dat
380c--------------------------------------------
381      if (g2d_premier) then
382        open (47,file='g2d.dat'
383clmd     &       ,form='unformatted',access='direct',recl=4)
384     &        ,form='unformatted',access='direct',recl=1
385     &        ,status='unknown')
386        g2d_irec=0
387        g2d_appel=0
388        g2d_premier=.false.
389      endif
390        g2d_appel = g2d_appel+1
391
392c--------------------------------------------
393c   Sortie g2d des xi proches + distants
394c--------------------------------------------
395cl                               if (nflev .NE. 500) then
396      do ja = 1,nuco2
397        do j = 0,nlaylte+1
398          do i = 0,nlaylte+1
399            g2d_irec=g2d_irec+1
400            reel4 = xi(ig1d,ja,i,j)
401            write(47,rec=g2d_irec) reel4
402          enddo
403        enddo
404      enddo
405cl                               endif
406
407c------------------------------------------------------
408c   Writeg2d des ksidb
409c------------------------------------------------------
410      do ja = 1,nuco2
411c       ja=1
412        do j = 0,nlaylte+1
413          do i = 0,nlaylte+1
414            g2d_irec=g2d_irec+1
415            reel4 = ksidb(ig,ja,i,j)
416            write(47,rec=g2d_irec) reel4
417          enddo
418        enddo
419      enddo
420
421      do j = 0,nlaylte+1
422        do i = 0,nlaylte+1
423          g2d_irec=g2d_irec+1
424          reel4 = ksidb(ig,3,i,j)
425          write(47,rec=g2d_irec) reel4
426        enddo
427      enddo
428
429c------------------------------------------------------
430c  Writeg2d dpsgcp
431c------------------------------------------------------
432
433        do j = 1 , nlaylte
434          do i = 0 , nlaylte+1
435            dpsgcp(i,j) = dp(ig,j) / gcp
436          enddo
437        enddo
438
439        do i = 0 , nlaylte+1
440c         dpsgcp(i,0) = 0.0002  ! (rapport ~ entre 1000 et 10000 pour le sol)
441          dpsgcp(i,0) = 1.      ! (pour regler l'echelle des sorties)
442          dpsgcp(i,nlaylte+1) = 0.
443        enddo
444
445c     print*
446c     print*,'gcp: ',gcp
447c     print*
448c       do i = 0 , nlaylte+1
449c     print*,i,'dp: ',dp(ig,i)
450c       enddo
451c     print*
452c       do i = 0 , nlaylte+1
453c     print*,i,'dpsgcp: ',dpsgcp(i,1)
454c       enddo
455 
456      do j = 0,nlaylte+1
457        do i = 0,nlaylte+1
458          g2d_irec=g2d_irec+1
459          reel4 = dpsgcp(i,j)
460          write(47,rec=g2d_irec) reel4
461        enddo
462      enddo
463
464c------------------------------------------------------
465c  Writeg2d temperature
466c------------------------------------------------------
467
468        do j = 1 , nlaylte
469          do i = 0 , nlaylte+1
470            temp(i,j) = tlay(ig,j)
471          enddo
472        enddo
473
474        do i = 0 , nlaylte+1
475          temp(i,0) = tlev(ig,1)+dt0(ig)     ! temperature surface
476          temp(i,nlaylte+1) = 0.               ! temperature espace  (=0)
477        enddo
478
479      do j = 0,nlaylte+1
480        do i = 0,nlaylte+1
481          g2d_irec=g2d_irec+1
482          reel4 = temp(i,j)
483          write(47,rec=g2d_irec) reel4
484        enddo
485      enddo
486
487        write(76,*) 'ig1d, ig, ig0', ig1d, ig, ig0
488        write(76,*) 'nlaylte', nlaylte
489        write(76,*) 'nflev', nflev
490        write(76,*) 'kdlon', kdlon
491        write(76,*) 'ndlo2', ndlo2
492        write(76,*) 'ndlon', ndlon
493      do ja=1,4
494        write(76,*) 'bsurf', ja, bsurf(ig,ja)
495        write(76,*) 'btop', ja, btop(ig,ja)
496
497        do j=1,nlaylte+1
498          write(76,*) 'blev', ja, j, blev(ig,ja,j)
499        enddo
500
501        do j=1,nlaylte
502          write(76,*) 'blay', ja, j, blay(ig,ja,j)
503        enddo
504
505        do j=1,2*nlaylte
506          write(76,*) 'dbsublay', ja, j, dbsublay(ig,ja,j)
507        enddo
508      enddo
509
510      endif
511c************************************************************************
512c#endif
513      endif  !   callg2d
514
515      return
516      end
Note: See TracBrowser for help on using the repository browser.