source: trunk/LMDZ.MARS/libf/phymars/lwxn.F @ 1009

Last change on this file since 1009 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: 12.4 KB
Line 
1      subroutine lwxn ( ig0,kdlon,kflev
2     .                , dp
3     .                , aer_t,co2_u,co2_up)
4
5c----------------------------------------------------------------------
6c     LWXN   computes transmission function and exchange coefficiants
7c                        for neighbours layers
8c                          (co2 / aerosols)
9c                       (bands 1 and 2 of co2)
10c----------------------------------------------------------------------
11
12c              |---|---|---|---|---|---|---|---|
13c   kflev+1    |   |   |   |   |   |   |   | 0 |  (space)
14c              |---|---|---|---|---|---|---|---|
15c    kflev     |   |   |   |   |   |***| 0 |   |
16c              |---|---|---|---|---|---|---|---|
17c     ...      |   |   |   |   |***| 0 |***|   |
18c              |---|---|---|---|---|---|---|---|
19c      4       |   |   |   |***| 0 |***|   |   |
20c              |---|---|---|---|---|---|---|---|
21c      3       |   |   |***| 0 |***|   |   |   |
22c              |---|---|---|---|---|---|---|---|
23c      2       |   |***| 0 |***|   |   |   |   |
24c              |---|---|---|---|---|---|---|---|
25c      1       |   | 0 |***|   |   |   |   |   |
26c              |---|---|---|---|---|---|---|---|
27c      0       | 0 |   |   |   |   |   |   |   |  (ground)
28c              |---|---|---|---|---|---|---|---|
29c                0   1   2   3   4  ...  k |k+1
30c             (ground)                    (space)
31
32c  (*)  xi computed in this subroutine
33c----------------------------------------------------------------------
34c                                                                       
35c  *********************************************************** nj=1
36c
37c
38c                                       sublayer    1                 
39c
40c
41c                                ----------------------------- nj=2
42c
43c                                       sublayer    2                 
44c
45c    - - - -  LAYER j - - - - -  ----------------------------- nj=3
46c
47c                                       sublayer    3                 
48c                                ----------------------------- nj=4
49c                                       sublayer  ncouche                 
50c  *********************************************************** ni=nj=ncouche+1
51c                                       sublayer  ncouche                 
52c                                ----------------------------- ni=4
53c                                       sublayer    3                 
54c
55c    - - - -  LAYER i - - - - -  ----------------------------- ni=3
56c
57c                                       sublayer    2                 
58c                     
59c                                ----------------------------- ni=2
60c
61c
62c                                       sublayer    1                 
63c
64c
65c  *********************************************************** ni=1
66c                                                                       
67c-----------------------------------------------------------------------
68c ATTENTION AUX UNITES:
69c le facteur 10*g fait passer des kg m-2 aux g cm-2
70c-----------------------------------------------------------------------
71
72      implicit none
73
74#include "dimensions.h"
75#include "dimphys.h"
76#include "dimradmars.h"
77 
78#include "yomlw.h"
79#include "callkeys.h"
80
81c----------------------------------------------------------------------
82c         0.1   arguments
83c               ---------
84c                                                            inputs:
85c                                                            -------
86      integer ig0
87      integer kdlon     ! part of ngrid
88      integer kflev     ! part of nalyer
89
90      real dp (ndlo2,kflev)              ! layer pressure thickness (Pa)
91
92      real aer_t (ndlo2,nuco2,kflev+1)   ! transmission (aer)
93      real co2_u (ndlo2,nuco2,kflev+1)   ! absorber amounts (co2)
94      real co2_up (ndlo2,nuco2,kflev+1)  ! idem scaled by the pressure (co2)
95
96c----------------------------------------------------------------------
97c         0.2   local arrays
98c               ------------
99
100      integer ja,jl,jk,nlmd,ni,nj
101
102      integer nmax
103      parameter (nmax=50)               ! max: 50 sublayers
104
105      real cn (nmax), cb (nmax)
106
107      real zu_layer_i (ndlon,nuco2)
108      real zup_layer_i (ndlon,nuco2)
109      real zt_aer_layer_i (ndlon,nuco2)
110
111      real zu_layer_j (ndlon,nuco2)
112      real zup_layer_j (ndlon,nuco2)
113      real zt_aer_layer_j (ndlon,nuco2)
114
115      real zu (ndlon,nuco2)
116      real zup (ndlon,nuco2)
117      real zt_co2 (ndlon,nuco2)
118      real zt_aer (ndlon,nuco2)
119
120      real zu_i (ndlon,nuco2,nmax+1)
121      real zup_i (ndlon,nuco2,nmax+1)
122      real zu_j (ndlon,nuco2,nmax+1)
123      real zup_j (ndlon,nuco2,nmax+1)
124      real zt_aer_i (ndlon,nuco2,nmax+1)
125      real zt_aer_j (ndlon,nuco2,nmax+1)
126
127      real trans (ndlon,nuco2,nmax+1,nmax+1)
128      real ksi (ndlon,nuco2,nflev-1)
129
130c----------------------------------------------------------------------
131c         0.3   Initialisation
132c               --------------
133
134      jk=ncouche+1
135      do ja = 1 ,nuco2
136        do jl = 1 , kdlon
137          zu_i (jl,ja,jk) = 0.
138          zup_i (jl,ja,jk) = 0.
139          zu_j (jl,ja,jk) = 0.
140          zup_j (jl,ja,jk) = 0.
141          zt_aer_i (jl,ja,jk) = 1.
142          zt_aer_j (jl,ja,jk) = 1.
143        enddo
144      enddo
145
146      if (linear) then
147
148      do nlmd = 1 ,ncouche
149        cn(nlmd)=(1.0/ncouche)
150        cb(nlmd)=(ncouche-nlmd+0.5)/ncouche
151c       print*,nlmd,cb(nlmd),cn(nlmd)
152      enddo
153
154      else
155
156      do nlmd = 1 ,ncouche-1
157        cn(nlmd)=(1-alphan)*alphan**(nlmd-1)
158        cb(nlmd)=0.5*(1+alphan)*alphan**(nlmd-1)
159      enddo
160      cn(ncouche)=alphan**(ncouche-1)
161      cb(ncouche)=0.5*alphan**(ncouche-1)
162
163      endif
164
165c test
166      if (nmax .LT. ncouche) then
167        print*,'!!!!! ATTENTION !!!!! '
168        print*,'probleme dans lwxn.F'
169        print*,' nmax=',nmax,'  < ncouche=',ncouche
170        call exit(1)
171      endif
172
173c----------------------------------------------------------------------
174      do jk = 1 , nlaylte-1
175c----------------------------------------------------------------------
176c         1.0   (co2) amount and (aer) transmission for all sublayers 
177c               ----------------------------------------------------
178
179        do ja = 1 , nuco2
180          do jl = 1 , kdlon
181
182c                                                        layer i (down)
183c                                                        -------------
184      zu_layer_i(jl,ja) =  co2_u(jl,ja,jk)  - co2_u(jl,ja,jk+1)
185      zup_layer_i(jl,ja) = co2_up(jl,ja,jk) - co2_up(jl,ja,jk+1)
186      zt_aer_layer_i(jl,ja) = aer_t(jl,ja,jk)
187     .                       / aer_t(jl,ja,jk+1)
188                                 
189      do nlmd=1,ncouche
190        zu_i(jl,ja,nlmd)=cn(nlmd)*zu_layer_i(jl,ja)
191        zup_i(jl,ja,nlmd)=cn(nlmd)*zup_layer_i(jl,ja)
192        zt_aer_i(jl,ja,nlmd)=zt_aer_layer_i(jl,ja)**cn(nlmd)
193      enddo
194
195c                                                          layer j (up)
196c                                                          ------------
197      zu_layer_j(jl,ja) =  co2_u(jl,ja,jk+1)  - co2_u(jl,ja,jk+2)
198      zup_layer_j(jl,ja) = co2_up(jl,ja,jk+1) - co2_up(jl,ja,jk+2)
199      zt_aer_layer_j(jl,ja) = aer_t(jl,ja,jk+1)
200     .                       / aer_t(jl,ja,jk+2)
201       
202      do nlmd=1,ncouche
203        zu_j(jl,ja,nlmd)=cn(nlmd)*zu_layer_j(jl,ja)
204        zup_j(jl,ja,nlmd)=cn(nlmd)*zup_layer_j(jl,ja)
205        zt_aer_j(jl,ja,nlmd)=zt_aer_layer_j(jl,ja)**cn(nlmd)
206      enddo
207
208          enddo
209        enddo
210
211c----------------------------------------------------------------------
212c         2.0   transmissions between all sublayers
213c               ------------------------------------
214
215        do ni = 1 ,ncouche+1
216
217            do ja = 1 ,nuco2
218              do jl = 1 , kdlon
219      zu(jl,ja)=0.   
220      zup(jl,ja)=0.   
221      zt_aer(jl,ja)=1.   
222                       
223      do nlmd=ni,ncouche+1
224        zu(jl,ja)=zu(jl,ja)+zu_i(jl,ja,nlmd)   
225        zup(jl,ja)=zup(jl,ja)+zup_i(jl,ja,nlmd)   
226        zt_aer(jl,ja)=zt_aer(jl,ja)*zt_aer_i(jl,ja,nlmd)   
227      enddo
228              enddo
229            enddo
230
231      call lwtt(kdlon,zu,zup,nuco2,zt_co2)
232                         
233        do ja = 1 ,nuco2
234          do jl = 1 , kdlon
235            trans(jl,ja,ni,ncouche+1)=zt_co2(jl,ja)*zt_aer(jl,ja)
236          enddo
237        enddo
238                   
239c on ajoute la couche J
240            do ja = 1 ,nuco2
241              do jl = 1 , kdlon
242      zu(jl,ja)=zu(jl,ja)+zu_layer_j(jl,ja)
243      zup(jl,ja)=zup(jl,ja)+zup_layer_j(jl,ja)
244      zt_aer(jl,ja)=zt_aer(jl,ja)*zt_aer_layer_j(jl,ja)
245              enddo
246            enddo
247                     
248      call lwtt(kdlon,zu,zup,nuco2,zt_co2)
249                         
250        do ja = 1 ,nuco2
251          do jl = 1 , kdlon
252            trans(jl,ja,ni,1)=zt_co2(jl,ja)*zt_aer(jl,ja)
253          enddo
254        enddo
255                   
256        enddo
257ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
258       
259        do nj = 1 ,ncouche+1
260
261            do ja = 1 ,nuco2
262              do jl = 1 , kdlon
263      zu(jl,ja)=0.   
264      zup(jl,ja)=0.   
265      zt_aer(jl,ja)=1.   
266                       
267      do nlmd=nj,ncouche+1
268        zu(jl,ja)=zu(jl,ja)+zu_j(jl,ja,nlmd)   
269        zup(jl,ja)=zup(jl,ja)+zup_j(jl,ja,nlmd)   
270        zt_aer(jl,ja)=zt_aer(jl,ja)*zt_aer_j(jl,ja,nlmd)   
271      enddo
272              enddo
273            enddo
274                     
275      call lwtt(kdlon,zu,zup,nuco2,zt_co2)
276                         
277        do ja = 1 ,nuco2
278          do jl = 1 , kdlon
279            trans(jl,ja,ncouche+1,nj)=zt_co2(jl,ja)*zt_aer(jl,ja)
280          enddo
281        enddo
282
283c on ajoute la couche I
284            do ja = 1 ,nuco2
285              do jl = 1 , kdlon
286      zu(jl,ja)=zu(jl,ja)+zu_layer_i(jl,ja)
287      zup(jl,ja)=zup(jl,ja)+zup_layer_i(jl,ja)
288      zt_aer(jl,ja)=zt_aer(jl,ja)*zt_aer_layer_i(jl,ja)
289              enddo
290            enddo
291                     
292      call lwtt(kdlon,zu,zup,nuco2,zt_co2)
293                         
294        do ja = 1 ,nuco2
295          do jl = 1 , kdlon
296            trans(jl,ja,1,nj)=zt_co2(jl,ja)*zt_aer(jl,ja)
297          enddo
298        enddo
299                   
300        enddo
301       
302c----------------------------------------------------------------------
303c         3.0   global exchange coefficiant between neigthbours
304c               -----------------------------------------------
305       
306        do ja = 1 ,nuco2
307          do jl = 1 , kdlon
308            ksi(jl,ja,jk) = 0.
309          enddo
310        enddo
311
312        do ni = 1 ,ncouche
313          do ja = 1 ,nuco2
314            do jl = 1 , kdlon
315
316      ksi(jl,ja,jk)=ksi(jl,ja,jk) +
317     .             ( trans(jl,ja,ni+1,ncouche+1)
318     .             - trans(jl,ja,ni,ncouche+1)
319     .             - trans(jl,ja,ni+1,1)
320     .             + trans(jl,ja,ni,1) )
321     .             * (cb(ni)*dp(jl,jk)) * 2   
322     .             /  (dp(jl,jk) + dp(jl,jk+1))       !!!!!!!!!!!!!!!!!!!
323
324            enddo
325          enddo
326        enddo
327
328        do nj = 1 ,ncouche
329          do ja = 1 ,nuco2
330            do jl = 1 , kdlon
331
332      ksi(jl,ja,jk)=ksi(jl,ja,jk) +
333     .             ( trans(jl,ja,ncouche+1,nj+1)
334     .             - trans(jl,ja,1,nj+1)
335     .             - trans(jl,ja,ncouche+1,nj)
336     .             + trans(jl,ja,1,nj) )
337     .             * (cb(nj)*dp(jl,jk+1)) * 2
338     .             /  (dp(jl,jk) + dp(jl,jk+1))       !!!!!!!!!!!!!!!!!!!
339
340            enddo
341          enddo
342        enddo
343
344        do ja = 1 ,nuco2
345          do jl = 1 , kdlon
346            xi(ig0+jl,ja,jk,jk+1) = ksi(jl,ja,jk)
347     .                            + xi_emis(ig0+jl,ja,jk)
348
349c                                                        ksi reciprocity
350c                                                        ---------------
351            xi(ig0+jl,ja,jk+1,jk) = xi(ig0+jl,ja,jk,jk+1)
352          enddo
353        enddo
354
355c----------------------------------------------------------------------
356c         4.0   Special treatment for ground
357c               ----------------------------
358
359
360      if (jk .EQ. 1) then
361
362        do ja = 1 ,nuco2
363          do jl = 1 , kdlon
364            xi_ground(ig0+jl,ja)=0.
365          enddo
366        enddo
367
368        do ni = 1 ,ncouche
369          do ja = 1 ,nuco2
370              do jl = 1 , kdlon
371
372      xi_ground(ig0+jl,ja) = xi_ground(ig0+jl,ja)
373     .                     + ( trans(jl,ja,ni+1,ncouche+1)
374     .                        -trans(jl,ja,ni,ncouche+1))
375     .                     * 2 * cb(ni)
376            enddo
377          enddo
378        enddo
379
380      endif
381
382c----------------------------------------------------------------------
383      enddo    !  boucle sur jk
384c----------------------------------------------------------------------
385      return
386      end
387
388
389
Note: See TracBrowser for help on using the repository browser.