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

Last change on this file since 2800 was 1917, checked in by emillour, 7 years ago

Mars GCM:
Code cleanup: get rid of routine "zerophys".
EM

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