source: trunk/LMDZ.MARS/libf/phymars/lwxd.F

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

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

File size: 9.2 KB
Line 
1      subroutine lwxd (ig0,kdlon,kflev,emis
2     .                ,aer_t,co2_u,co2_up)
3
4c----------------------------------------------------------------------
5c     LWXD   computes transmission function and exchange coefficiants
6c                        for distant layers
7c                          (co2 / aerosols)
8c                       (bands 1 and 2 of co2)
9c----------------------------------------------------------------------
10
11c              |---|---|---|---|---|---|---|---|
12c   kflev+1    |   |   |   |   |   |   |   | 0 |  (space)
13c              |---|---|---|---|---|---|---|---|
14c    kflev     |   |***|***|***|***|   | 0 |   |
15c              |---|---|---|---|---|---|---|---|
16c     ...      |   |***|***|***|   | 0 |   |   |
17c              |---|---|---|---|---|---|---|---|
18c      4       |   |***|***|   | 0 |   |***|   |
19c              |---|---|---|---|---|---|---|---|
20c      3       |   |***|   | 0 |   |***|***|   |
21c              |---|---|---|---|---|---|---|---|
22c      2       |   |   | 0 |   |   |***|***|   |
23c              |---|---|---|---|---|---|---|---|
24c      1       |   | 0 |   |   |***|***|***|   |
25c              |---|---|---|---|---|---|---|---|
26c      0       | 0 |   |   |***|***|***|***|   |  (ground)
27c              |---|---|---|---|---|---|---|---|
28c                0   1   2   3   4  ...  k |k+1
29c             (ground)                    (space)
30
31c  (*)  xi computed in this subroutine
32c----------------------------------------------------------------------
33
34      use dimradmars_mod, only: ndlon, nuco2, nflev, ndlo2
35      use yomlw_h, only: nlaylte, xi, xi_emis
36      implicit none
37
38      include "callkeys.h"
39
40c----------------------------------------------------------------------
41c         0.1   arguments
42c               ---------
43c                                                            inputs:
44c                                                            -------
45      integer,intent(in) :: ig0
46      integer,intent(in) :: kdlon      ! part of ngrid
47      integer,intent(in) :: kflev      ! part of nalyer
48 
49      real,intent(in) :: emis (ndlo2) ! surface emissivity
50      real,intent(in) :: aer_t (ndlo2,nuco2,kflev+1) ! transmission (aer)
51      real,intent(in) :: co2_u (ndlo2,nuco2,kflev+1) ! absorber amounts (co2)
52      real,intent(in) :: co2_up (ndlo2,nuco2,kflev+1) ! idem scaled by the pressure (co2)
53
54c----------------------------------------------------------------------
55c         0.2   local arrays
56c               ------------
57
58      integer ja,jl,jk,jkk
59      real zu (ndlon,nuco2)
60      real zup (ndlon,nuco2)
61      real zt_co2 (ndlon,nuco2)
62      real zt_aer (ndlon,nuco2)
63
64      real ksi (ndlon,nuco2,0:nflev+1,0:nflev+1)
65      real ksi_emis (ndlon,nuco2,0:nflev+1,0:nflev+1)
66      real trans (ndlon,nuco2,0:nflev+1,0:nflev+1)
67      real trans_emis (ndlon,nuco2,0:nflev+1,0:nflev+1)
68
69c----------------------------------------------------------------------
70      ksi_emis(:,:,:,:)=0
71c----------------------------------------------------------------------
72c         1.0   Transmission functions
73c               ----------------------
74
75c----------------------------------------------------------------------
76c        1.1   Direct transmission
77c              -------------------
78
79      do jk = 1 , nlaylte+1
80        do jkk = jk , nlaylte+1
81
82          do ja = 1 , nuco2
83            do jl = 1 , kdlon
84c                                                                   co2
85c                                                                   ---
86              zu(jl,ja) =  co2_u(jl,ja,jk)  - co2_u(jl,ja,jkk)
87              zup(jl,ja) = co2_up(jl,ja,jk) - co2_up(jl,ja,jkk)
88c                                                                   aer
89c                                                                   ---
90              zt_aer(jl,ja)= aer_t(jl,ja,jk)
91     .                      /aer_t(jl,ja,jkk)
92
93            enddo
94          enddo
95
96          call lwtt(kdlon,zu,zup,nuco2,zt_co2)
97c                                                            co2 and aer
98c                                                            -----------
99          do ja = 1 , nuco2
100            do jl = 1 , kdlon
101              trans(jl,ja,jk,jkk) = zt_co2(jl,ja) * zt_aer(jl,ja)
102            enddo
103          enddo
104c                                                       trans reciprocity
105c                                                       -----------------
106          do ja = 1 , nuco2
107            do jl = 1 , kdlon
108              trans(jl,ja,jkk,jk) = trans(jl,ja,jk,jkk)
109c           if (trans(jl,ja,jk,jkk) .LT. 0 ) then
110c           print*,'trans bande',ja,jk,jkk,trans(jl,ja,jk,jkk)
111c           endif
112c           if (trans(jl,ja,jk,jkk) .GT. 1) then
113c           print*,'trans bande',ja,jk,jkk,trans(jl,ja,jk,jkk)
114c           trans(jl,ja,jk,jkk)=1
115c           print*,'trans bande',ja,jk,jkk,trans(jl,ja,jk,jkk)
116c           endif
117
118            enddo
119          enddo
120     
121        enddo
122      enddo
123
124c----------------------------------------------------------------------
125c         1.2   Transmission with reflexion
126c               ---------------------------
127
128      do jk = 1 , nlaylte+1
129        do jkk = jk , nlaylte+1
130
131      if (callemis) then
132          do ja = 1 , nuco2
133            do jl = 1 , kdlon
134c                                                                   co2
135c                                                                   ---
136              zu(jl,ja)  = 2 * co2_u(jl,ja,1)  - co2_u(jl,ja,jk)
137     .                                         - co2_u(jl,ja,jkk)
138              zup(jl,ja) = 2 * co2_up(jl,ja,1) - co2_up(jl,ja,jk)
139     .                                         - co2_up(jl,ja,jkk)
140c                                                                   aer
141c                                                                   ---
142                zt_aer(jl,ja)  = aer_t(jl,ja,1)
143     .                         * aer_t(jl,ja,1)
144     .                         / aer_t(jl,ja,jk)
145     .                         / aer_t(jl,ja,jkk)
146            enddo
147          enddo
148
149          call lwtt(kdlon,zu,zup,nuco2,zt_co2)
150c                                                            co2 and aer
151c                                                            -----------
152          do ja = 1 , nuco2
153            do jl = 1 , kdlon
154              trans_emis(jl,ja,jk,jkk) = zt_co2(jl,ja)
155     .                                 * zt_aer(jl,ja)
156            enddo
157          enddo
158
159      else
160
161          do ja = 1 , nuco2
162            do jl = 1 , kdlon
163              trans_emis(jl,ja,jk,jkk) = 1.
164            enddo
165          enddo
166
167      endif
168c                                                       trans reciprocity
169c                                                       -----------------
170          do ja = 1 , nuco2
171            do jl = 1 , kdlon
172              trans_emis(jl,ja,jkk,jk) = trans_emis(jl,ja,jk,jkk)
173c         if (trans_emis(jl,ja,jk,jkk) .LT. 0
174c    .                .OR.  trans_emis(jl,ja,jk,jkk) .GT. 1) then
175c     print*,'trans_emis bande',ja,jk,jkk,trans_emis(jl,ja,jk,jkk)
176c         endif
177            enddo
178          enddo
179
180        enddo
181      enddo
182
183c----------------------------------------------------------------------
184c         2.0   Exchange Coefficiants
185c               ---------------------
186
187      do jk = 1 , nlaylte-2
188        do jkk = jk+2 , nlaylte
189          do ja = 1 , nuco2
190            do jl = 1 , kdlon
191
192      ksi(jl,ja,jk,jkk) =
193     .            trans(jl,ja,jk+1,jkk)   - trans(jl,ja,jk,jkk)
194     .          - trans(jl,ja,jk+1,jkk+1) + trans(jl,ja,jk,jkk+1)
195
196      ksi_emis(jl,ja,jk,jkk) =
197     .   trans_emis(jl,ja,jk,jkk)   - trans_emis(jl,ja,jk+1,jkk)
198     . - trans_emis(jl,ja,jk,jkk+1) + trans_emis(jl,ja,jk+1,jkk+1)
199
200c       if (ksi(jl,ja,jk,jkk) .LT. 0 ) then
201c           print*,'ksi bande',ja,jk,jkk,ksi(jl,ja,jk,jkk)
202c           ksi(jl,ja,jk,jkk)=0
203c           print*,'ksi bande',ja,jk,jkk,ksi(jl,ja,jk,jkk)
204c       endif
205c       if (ksi(jl,ja,jk,jkk) .GT. 1) then
206c           print*,'ksi bande',ja,jk,jkk,ksi(jl,ja,jk,jkk)
207c           ksi(jl,ja,jk,jkk)=1
208c           print*,'ksi bande',ja,jk,jkk,ksi(jl,ja,jk,jkk)
209c       endif
210
211c       if (ksi_emis(jl,ja,jk,jkk) .LT. 0
212c    .                .OR.  ksi_emis(jl,ja,jk,jkk) .GT. 1) then
213c     print*,'ksi_emis bande',ja,jk,jkk,ksi_emis(jl,ja,jk,jkk)
214c       endif
215
216      xi(ig0+jl,ja,jk,jkk) = ksi(jl,ja,jk,jkk)
217     .      + ksi_emis(jl,ja,jk,jkk) * (1 - emis(jl))
218
219c                                                        ksi reciprocity
220c                                                        ---------------
221      ksi(jl,ja,jkk,jk)      = ksi(jl,ja,jk,jkk)
222      ksi_emis(jl,ja,jkk,jk) = ksi_emis(jl,ja,jk,jkk)
223      xi(ig0+jl,ja,jkk,jk)   = xi(ig0+jl,ja,jk,jkk)
224
225            enddo
226          enddo
227        enddo
228      enddo
229
230c----------------------------------------------------------------------
231c         2.1   Save xi_emis for neighbours (lwxn.F)
232c               -----------------------------------
233
234      do jk = 1 , nlaylte-1
235        do ja = 1 , nuco2
236          do jl = 1 , kdlon
237
238c     ksi_emis(jl,ja,jk,jk+1) =
239c    .   trans_emis(jl,ja,jk,jk+1)   - trans_emis(jl,ja,jk+1,jk+1)
240c    . - trans_emis(jl,ja,jk,jk+2) + trans_emis(jl,ja,jk+1,jk+2)
241
242      xi_emis(ig0+jl,ja,jk) =
243     .                 ksi_emis(jl,ja,jk,jk+1) * (1-emis(jl))
244
245          enddo
246        enddo
247      enddo
248
249c----------------------------------------------------------------------
250      end
Note: See TracBrowser for help on using the repository browser.