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

Last change on this file since 3118 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

  • Remove all reference/use of nlayermx and dimphys.h
  • Made use of automatic arrays whenever arrays are needed with dimension nlayer
  • Remove lots of obsolete reference to dimensions.h
  • Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

  • Deleted param.h and put contents into module param_v4_h
  • Adapted testphys1d, newstart, etc...
  • Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

  • Optimized computations in paramfoto_compact (twice less dlog10 calculations)
  • Checked consistency before/after modification in debug mode
  • Checked performance is not impacted (same as before)
File size: 12.5 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      use dimradmars_mod, only: ndlo2, nuco2, ndlon, nflev
73      use yomlw_h, only: nlaylte, xi, xi_ground, xi_emis
74      implicit none
75
76#include "callkeys.h"
77
78c----------------------------------------------------------------------
79c         0.1   arguments
80c               ---------
81c                                                            inputs:
82c                                                            -------
83      integer ig0
84      integer kdlon     ! part of ngrid
85      integer kflev     ! part of nalyer
86
87      real dp (ndlo2,kflev)              ! layer pressure thickness (Pa)
88
89      real aer_t (ndlo2,nuco2,kflev+1)   ! transmission (aer)
90      real co2_u (ndlo2,nuco2,kflev+1)   ! absorber amounts (co2)
91      real co2_up (ndlo2,nuco2,kflev+1)  ! idem scaled by the pressure (co2)
92
93c----------------------------------------------------------------------
94c         0.2   local arrays
95c               ------------
96
97      integer ja,jl,jk,nlmd,ni,nj
98
99      integer nmax
100      parameter (nmax=50)               ! max: 50 sublayers
101
102      real cn (nmax), cb (nmax)
103
104      real zu_layer_i (ndlon,nuco2)
105      real zup_layer_i (ndlon,nuco2)
106      real zt_aer_layer_i (ndlon,nuco2)
107
108      real zu_layer_j (ndlon,nuco2)
109      real zup_layer_j (ndlon,nuco2)
110      real zt_aer_layer_j (ndlon,nuco2)
111
112      real zu (ndlon,nuco2)
113      real zup (ndlon,nuco2)
114      real zt_co2 (ndlon,nuco2)
115      real zt_aer (ndlon,nuco2)
116
117      real zu_i (ndlon,nuco2,nmax+1)
118      real zup_i (ndlon,nuco2,nmax+1)
119      real zu_j (ndlon,nuco2,nmax+1)
120      real zup_j (ndlon,nuco2,nmax+1)
121      real zt_aer_i (ndlon,nuco2,nmax+1)
122      real zt_aer_j (ndlon,nuco2,nmax+1)
123
124      real trans (ndlon,nuco2,nmax+1,nmax+1)
125      real ksi (ndlon,nuco2,nflev-1)
126
127c----------------------------------------------------------------------
128c         0.3   Initialisation
129c               --------------
130
131      jk=ncouche+1
132      do ja = 1 ,nuco2
133        do jl = 1 , kdlon
134          zu_i (jl,ja,jk) = 0.
135          zup_i (jl,ja,jk) = 0.
136          zu_j (jl,ja,jk) = 0.
137          zup_j (jl,ja,jk) = 0.
138          zt_aer_i (jl,ja,jk) = 1.
139          zt_aer_j (jl,ja,jk) = 1.
140        enddo
141      enddo
142
143      if (linear) then
144
145      do nlmd = 1 ,ncouche
146        cn(nlmd)=(1.0/ncouche)
147        cb(nlmd)=(ncouche-nlmd+0.5)/ncouche
148c       print*,nlmd,cb(nlmd),cn(nlmd)
149      enddo
150
151      else
152
153      do nlmd = 1 ,ncouche-1
154        cn(nlmd)=(1-alphan)*alphan**(nlmd-1)
155        cb(nlmd)=0.5*(1+alphan)*alphan**(nlmd-1)
156      enddo
157      cn(ncouche)=alphan**(ncouche-1)
158      cb(ncouche)=0.5*alphan**(ncouche-1)
159
160      endif
161
162c test
163      if (nmax .LT. ncouche) then
164        print*,'!!!!! ATTENTION !!!!! '
165        print*,'probleme dans lwxn.F'
166        print*,' nmax=',nmax,'  < ncouche=',ncouche
167        call exit(1)
168      endif
169
170c----------------------------------------------------------------------
171      do jk = 1 , nlaylte-1
172c----------------------------------------------------------------------
173c         1.0   (co2) amount and (aer) transmission for all sublayers 
174c               ----------------------------------------------------
175
176        do ja = 1 , nuco2
177          do jl = 1 , kdlon
178
179c                                                        layer i (down)
180c                                                        -------------
181      zu_layer_i(jl,ja) =  co2_u(jl,ja,jk)  - co2_u(jl,ja,jk+1)
182      zup_layer_i(jl,ja) = co2_up(jl,ja,jk) - co2_up(jl,ja,jk+1)
183      zt_aer_layer_i(jl,ja) = aer_t(jl,ja,jk)
184     .                       / aer_t(jl,ja,jk+1)
185                                 
186      do nlmd=1,ncouche
187        zu_i(jl,ja,nlmd)=cn(nlmd)*zu_layer_i(jl,ja)
188        zup_i(jl,ja,nlmd)=cn(nlmd)*zup_layer_i(jl,ja)
189        zt_aer_i(jl,ja,nlmd)=zt_aer_layer_i(jl,ja)**cn(nlmd)
190      enddo
191
192c                                                          layer j (up)
193c                                                          ------------
194      zu_layer_j(jl,ja) =  co2_u(jl,ja,jk+1)  - co2_u(jl,ja,jk+2)
195      zup_layer_j(jl,ja) = co2_up(jl,ja,jk+1) - co2_up(jl,ja,jk+2)
196      zt_aer_layer_j(jl,ja) = aer_t(jl,ja,jk+1)
197     .                       / aer_t(jl,ja,jk+2)
198       
199      do nlmd=1,ncouche
200        zu_j(jl,ja,nlmd)=cn(nlmd)*zu_layer_j(jl,ja)
201        zup_j(jl,ja,nlmd)=cn(nlmd)*zup_layer_j(jl,ja)
202        zt_aer_j(jl,ja,nlmd)=zt_aer_layer_j(jl,ja)**cn(nlmd)
203      enddo
204
205          enddo
206        enddo
207
208c----------------------------------------------------------------------
209c         2.0   transmissions between all sublayers
210c               ------------------------------------
211
212        do ni = 1 ,ncouche+1
213
214            do ja = 1 ,nuco2
215              do jl = 1 , kdlon
216      zu(jl,ja)=0.   
217      zup(jl,ja)=0.   
218      zt_aer(jl,ja)=1.   
219                       
220      do nlmd=ni,ncouche+1
221        zu(jl,ja)=zu(jl,ja)+zu_i(jl,ja,nlmd)   
222        zup(jl,ja)=zup(jl,ja)+zup_i(jl,ja,nlmd)   
223        zt_aer(jl,ja)=zt_aer(jl,ja)*zt_aer_i(jl,ja,nlmd)   
224      enddo
225              enddo
226            enddo
227
228      call lwtt(kdlon,zu,zup,nuco2,zt_co2)
229                         
230        do ja = 1 ,nuco2
231          do jl = 1 , kdlon
232            trans(jl,ja,ni,ncouche+1)=zt_co2(jl,ja)*zt_aer(jl,ja)
233          enddo
234        enddo
235                   
236c on ajoute la couche J
237            do ja = 1 ,nuco2
238              do jl = 1 , kdlon
239      zu(jl,ja)=zu(jl,ja)+zu_layer_j(jl,ja)
240      zup(jl,ja)=zup(jl,ja)+zup_layer_j(jl,ja)
241      zt_aer(jl,ja)=zt_aer(jl,ja)*zt_aer_layer_j(jl,ja)
242              enddo
243            enddo
244                     
245      call lwtt(kdlon,zu,zup,nuco2,zt_co2)
246                         
247        do ja = 1 ,nuco2
248          do jl = 1 , kdlon
249            trans(jl,ja,ni,1)=zt_co2(jl,ja)*zt_aer(jl,ja)
250          enddo
251        enddo
252                   
253        enddo
254ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
255       
256        do nj = 1 ,ncouche+1
257
258            do ja = 1 ,nuco2
259              do jl = 1 , kdlon
260      zu(jl,ja)=0.   
261      zup(jl,ja)=0.   
262      zt_aer(jl,ja)=1.   
263                       
264      do nlmd=nj,ncouche+1
265        zu(jl,ja)=zu(jl,ja)+zu_j(jl,ja,nlmd)   
266        zup(jl,ja)=zup(jl,ja)+zup_j(jl,ja,nlmd)   
267        zt_aer(jl,ja)=zt_aer(jl,ja)*zt_aer_j(jl,ja,nlmd)   
268      enddo
269              enddo
270            enddo
271                     
272      call lwtt(kdlon,zu,zup,nuco2,zt_co2)
273                         
274        do ja = 1 ,nuco2
275          do jl = 1 , kdlon
276            trans(jl,ja,ncouche+1,nj)=zt_co2(jl,ja)*zt_aer(jl,ja)
277          enddo
278        enddo
279
280c on ajoute la couche I
281            do ja = 1 ,nuco2
282              do jl = 1 , kdlon
283      zu(jl,ja)=zu(jl,ja)+zu_layer_i(jl,ja)
284      zup(jl,ja)=zup(jl,ja)+zup_layer_i(jl,ja)
285      zt_aer(jl,ja)=zt_aer(jl,ja)*zt_aer_layer_i(jl,ja)
286              enddo
287            enddo
288                     
289      call lwtt(kdlon,zu,zup,nuco2,zt_co2)
290                         
291        do ja = 1 ,nuco2
292          do jl = 1 , kdlon
293            trans(jl,ja,1,nj)=zt_co2(jl,ja)*zt_aer(jl,ja)
294          enddo
295        enddo
296                   
297        enddo
298       
299c----------------------------------------------------------------------
300c         3.0   global exchange coefficiant between neigthbours
301c               -----------------------------------------------
302       
303        do ja = 1 ,nuco2
304          do jl = 1 , kdlon
305            ksi(jl,ja,jk) = 0.
306          enddo
307        enddo
308
309        do ni = 1 ,ncouche
310          do ja = 1 ,nuco2
311            do jl = 1 , kdlon
312
313      ksi(jl,ja,jk)=ksi(jl,ja,jk) +
314     .             ( trans(jl,ja,ni+1,ncouche+1)
315     .             - trans(jl,ja,ni,ncouche+1)
316     .             - trans(jl,ja,ni+1,1)
317     .             + trans(jl,ja,ni,1) )
318     .             * (cb(ni)*dp(jl,jk)) * 2   
319     .             /  (dp(jl,jk) + dp(jl,jk+1))       !!!!!!!!!!!!!!!!!!!
320
321            enddo
322          enddo
323        enddo
324
325        do nj = 1 ,ncouche
326          do ja = 1 ,nuco2
327            do jl = 1 , kdlon
328
329      ksi(jl,ja,jk)=ksi(jl,ja,jk) +
330     .             ( trans(jl,ja,ncouche+1,nj+1)
331     .             - trans(jl,ja,1,nj+1)
332     .             - trans(jl,ja,ncouche+1,nj)
333     .             + trans(jl,ja,1,nj) )
334     .             * (cb(nj)*dp(jl,jk+1)) * 2
335     .             /  (dp(jl,jk) + dp(jl,jk+1))       !!!!!!!!!!!!!!!!!!!
336
337            enddo
338          enddo
339        enddo
340
341        do ja = 1 ,nuco2
342          do jl = 1 , kdlon
343            xi(ig0+jl,ja,jk,jk+1) = ksi(jl,ja,jk)
344     .                            + xi_emis(ig0+jl,ja,jk)
345
346c                                                        ksi reciprocity
347c                                                        ---------------
348            xi(ig0+jl,ja,jk+1,jk) = xi(ig0+jl,ja,jk,jk+1)
349          enddo
350        enddo
351
352c----------------------------------------------------------------------
353c         4.0   Special treatment for ground
354c               ----------------------------
355
356
357      if (jk .EQ. 1) then
358
359        do ja = 1 ,nuco2
360          do jl = 1 , kdlon
361            xi_ground(ig0+jl,ja)=0.
362          enddo
363        enddo
364
365        do ni = 1 ,ncouche
366          do ja = 1 ,nuco2
367              do jl = 1 , kdlon
368
369      xi_ground(ig0+jl,ja) = xi_ground(ig0+jl,ja)
370     .                     + ( trans(jl,ja,ni+1,ncouche+1)
371     .                        -trans(jl,ja,ni,ncouche+1))
372     .                     * 2 * cb(ni)
373            enddo
374          enddo
375        enddo
376
377      endif
378
379c----------------------------------------------------------------------
380      enddo    !  boucle sur jk
381c----------------------------------------------------------------------
382      return
383      end
384
385
386
Note: See TracBrowser for help on using the repository browser.