source: trunk/mesoscale/LMDZ.MARS/libf_gcm/phymars/lwflux.F @ 57

Last change on this file since 57 was 57, checked in by aslmd, 14 years ago

mineur LMD_MM_MARS: ajout du GCM ancienne physique, systeme maintenant complet sur SVN (ne manque que la base de donnees d'etats initiaux)

File size: 15.3 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(ig,l)   ;   downward : flw_dn(ig,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
362c ig:   point pour les variables decoupees
363
364c#ifdef undim
365      if (callg2d) then
366
367      ig1d = ngridmx/2 + 1
368c     ig1d = ngridmx
369
370      if ((ig0+1).LE.ig1d .and. ig1d.LE.(ig0+kdlon)
371     .    .OR.  ngridmx.EQ.1   ) then
372
373          ig = ig1d-ig0
374        print*, 'Sortie g2d: ig1d, ig, ig0', ig1d, ig, ig0
375
376c--------------------------------------------
377c   Ouverture de g2d.dat
378c--------------------------------------------
379      if (g2d_premier) then
380        open (47,file='g2d.dat'
381clmd     &       ,form='unformatted',access='direct',recl=4)
382     &        ,form='unformatted',access='direct',recl=1
383     &        ,status='unknown')
384        g2d_irec=0
385        g2d_appel=0
386        g2d_premier=.false.
387      endif
388        g2d_appel = g2d_appel+1
389
390c--------------------------------------------
391c   Sortie g2d des xi proches + distants
392c--------------------------------------------
393cl                               if (nflev .NE. 500) then
394      do ja = 1,nuco2
395        do j = 0,nlaylte+1
396          do i = 0,nlaylte+1
397            g2d_irec=g2d_irec+1
398            reel4 = xi(ig1d,ja,i,j)
399            write(47,rec=g2d_irec) reel4
400          enddo
401        enddo
402      enddo
403cl                               endif
404
405c------------------------------------------------------
406c   Writeg2d des ksidb
407c------------------------------------------------------
408      do ja = 1,nuco2
409c       ja=1
410        do j = 0,nlaylte+1
411          do i = 0,nlaylte+1
412            g2d_irec=g2d_irec+1
413            reel4 = ksidb(ig,ja,i,j)
414            write(47,rec=g2d_irec) reel4
415          enddo
416        enddo
417      enddo
418
419      do j = 0,nlaylte+1
420        do i = 0,nlaylte+1
421          g2d_irec=g2d_irec+1
422          reel4 = ksidb(ig,3,i,j)
423          write(47,rec=g2d_irec) reel4
424        enddo
425      enddo
426
427c------------------------------------------------------
428c  Writeg2d dpsgcp
429c------------------------------------------------------
430
431        do j = 1 , nlaylte
432          do i = 0 , nlaylte+1
433            dpsgcp(i,j) = dp(ig,j) / gcp
434          enddo
435        enddo
436
437        do i = 0 , nlaylte+1
438c         dpsgcp(i,0) = 0.0002  ! (rapport ~ entre 1000 et 10000 pour le sol)
439          dpsgcp(i,0) = 1.      ! (pour regler l'echelle des sorties)
440          dpsgcp(i,nlaylte+1) = 0.
441        enddo
442
443c     print*
444c     print*,'gcp: ',gcp
445c     print*
446c       do i = 0 , nlaylte+1
447c     print*,i,'dp: ',dp(ig,i)
448c       enddo
449c     print*
450c       do i = 0 , nlaylte+1
451c     print*,i,'dpsgcp: ',dpsgcp(i,1)
452c       enddo
453 
454      do j = 0,nlaylte+1
455        do i = 0,nlaylte+1
456          g2d_irec=g2d_irec+1
457          reel4 = dpsgcp(i,j)
458          write(47,rec=g2d_irec) reel4
459        enddo
460      enddo
461
462c------------------------------------------------------
463c  Writeg2d temperature
464c------------------------------------------------------
465
466        do j = 1 , nlaylte
467          do i = 0 , nlaylte+1
468            temp(i,j) = tlay(ig,j)
469          enddo
470        enddo
471
472        do i = 0 , nlaylte+1
473          temp(i,0) = tlev(ig,1)+dt0(ig)     ! temperature surface
474          temp(i,nlaylte+1) = 0.               ! temperature espace  (=0)
475        enddo
476
477      do j = 0,nlaylte+1
478        do i = 0,nlaylte+1
479          g2d_irec=g2d_irec+1
480          reel4 = temp(i,j)
481          write(47,rec=g2d_irec) reel4
482        enddo
483      enddo
484
485        write(76,*) 'ig1d, ig, ig0', ig1d, ig, ig0
486        write(76,*) 'nlaylte', nlaylte
487        write(76,*) 'nflev', nflev
488        write(76,*) 'kdlon', kdlon
489        write(76,*) 'ndlo2', ndlo2
490        write(76,*) 'ndlon', ndlon
491      do ja=1,4
492        write(76,*) 'bsurf', ja, bsurf(ig,ja)
493        write(76,*) 'btop', ja, btop(ig,ja)
494
495        do j=1,nlaylte+1
496          write(76,*) 'blev', ja, j, blev(ig,ja,j)
497        enddo
498
499        do j=1,nlaylte
500          write(76,*) 'blay', ja, j, blay(ig,ja,j)
501        enddo
502
503        do j=1,2*nlaylte
504          write(76,*) 'dbsublay', ja, j, dbsublay(ig,ja,j)
505        enddo
506      enddo
507
508      endif
509c************************************************************************
510c#endif
511      endif  !   callg2d
512
513      return
514      end
Note: See TracBrowser for help on using the repository browser.