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

Last change on this file since 706 was 38, checked in by emillour, 14 years ago

Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs.
EM

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