source: trunk/LMDZ.MARS/libf/phymars/lwi.F @ 1465

Last change on this file since 1465 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: 5.2 KB
RevLine 
[38]1      subroutine lwi (ig0,kdlon,kflev
2     .                ,psi,zdblay,pdp
3     .                ,newpcolc )
4
5
[1047]6      use dimradmars_mod, only: ndlo2, ndlon, nflev, nir
7      use yomlw_h, only: gcp, nlaylte, xi
[1226]8      USE comcstfi_h
[38]9      implicit none
10
11#include "comg1d.h"
12#include "callkeys.h"
13 
14CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
15C
16C                             -   lwi    -   
17C
18C     PURPOSE:       Shema semi - implicite
19C
20C
21CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
22
23
24c************************************************************************
25c
26c        0.    Declarations
27c              ------------
28c
29c-------------------------------------------------------------------------
30c        0.1   Arguments
31c              ---------
32c
33 
34      integer ig0,kdlon,kflev
35
36      real    psi(ndlo2,kflev)
37     .     ,  zdblay(ndlo2,nir,kflev)
38     .     ,  pdp(ndlo2,kflev)
39
40
41      real    newpcolc(ndlo2,kflev)
42
43c-------------------------------------------------------------------------
44c        0.2   local arrays
45c              ------------
46c
47      real    di(ndlon,nflev)
48     .      , hi(ndlon,nflev)
49     .      , bi(ndlon,nflev)
50
51      real    ci(ndlon,nflev)
52     .      , ai(ndlon,nflev)
53      real    deltat
54
55      real   semit, denom
56
57      integer i, jl
58
59c************************************************************************
60c
61c        1.    Initialisations
62c              ---------------
63c
64c-----------------------------------------------------------------------
65 
66        deltat = dtphys * iradia
67c       print*,'SEMI = ',semi, '(expl:0  semi-implicite:0.5  impl:1)'
68        semit = semi * deltat
69c       semi = 0.
70
71c       print*,'dtphys,iradia,deltat,semit:',dtphys,iradia,deltat,semit
72c       print*,'g,cpp',g,cpp
73
74
75c************************************************************************
76c
77c        2.   
78c              ---------------
79c
80c-------------------------------------------------------------------------
81c        2.1   Calcul des di
82c              -------------
83c
84
85
86      do i = 1 , nlaylte-1
87        do jl = 1 , kdlon
88c     -------------------
89      di(jl,i) =  1 + semit * (g / pdp(jl,i) / cpp) * (
90     .    ( xi(ig0+jl,1,i,nlaylte+1)
91     .    + xi(ig0+jl,1,i,i+1)
92     .    + xi(ig0+jl,1,i,i-1) )
93     .    *    zdblay(jl,1,i)
94     .  + ( xi(ig0+jl,2,i,nlaylte+1)
95     .    + xi(ig0+jl,2,i,i+1)
96     .    + xi(ig0+jl,2,i,i-1) )
97     .    *    zdblay(jl,2,i)
98     .     )
99c     -------------------
100        enddo
101      enddo
102
103c couche nlaylte
104c ------------
105c      , on enleve i,i+1 sinon on a 2 fois le cooling2space
106
107      do jl = 1 , kdlon
108c     -------------------
109      di(jl,nlaylte) =  1 + semit * (g / pdp(jl,nlaylte) / cpp) * (
110     .    ( xi(ig0+jl,1,nlaylte,nlaylte+1)
111     .    + xi(ig0+jl,1,nlaylte,nlaylte-1) )
112     .    *    zdblay(jl,1,nlaylte)
113     .  + ( xi(ig0+jl,2,nlaylte,nlaylte+1)
114     .    + xi(ig0+jl,2,nlaylte,nlaylte-1) )
115     .    *    zdblay(jl,2,nlaylte)
116     .     )
117c     -------------------
118      enddo
119
120c-------------------------------------------------------------------------
121c        2.2   Calcul des hi
122c              -------------
123c
124
125      do i = 1 , nlaylte-1
126        do jl = 1 , kdlon
127c     -------------------
128      hi(jl,i) =    - semit * (g / pdp(jl,i) / cpp) *
129     .            (    xi(ig0+jl,1,i,i+1) * zdblay(jl,1,i+1)   
130     .               + xi(ig0+jl,2,i,i+1) * zdblay(jl,2,i+1)   )
131c     -------------------
132        enddo
133      enddo
134
135c-------------------------------------------------------------------------
136c        2.3   Calcul des bi
137c              -------------
138c
139
140
141      do i = 2 , nlaylte
142        do jl = 1 , kdlon
143c     -------------------
144      bi(jl,i) =   - semit * (g / pdp(jl,i) / cpp) *
145     .           (     xi(ig0+jl,1,i,i-1) * zdblay(jl,1,i-1)   
146     .               + xi(ig0+jl,2,i,i-1) * zdblay(jl,2,i-1)   )
147c     -------------------
148        enddo
149      enddo
150
151
152c couche 1
153c --------
154c  tant qu'on a pas un calcul propre de zdblay(0) qui tienne compte de
155c    la discontinuite de temperature au sol , on met  b1 = 0
156
157
158      do jl = 1 , kdlon
159        bi(jl,1) = 0
160      enddo
161
162c-------------------------------------------------------------------------
163c        2.4   
164c              -------------
165c
166
167c couche nlaylte
168c ------------
169
170      do jl = 1 , kdlon
171c     -------------------
172      ci(jl,nlaylte) = (gcp * psi(jl,nlaylte) / pdp(jl,nlaylte))
173     .                   / di(jl,nlaylte)
174
175      ai(jl,nlaylte) = - bi(jl,nlaylte) / di(jl,nlaylte)
176c     -------------------
177      enddo
178
179
180
181      do i = nlaylte-1 , 1 , -1
182        do jl = 1 , kdlon
183c     -------------------
184      denom = di(jl,i) + hi(jl,i) * ai(jl,i+1)
185
186      ci(jl,i) = (  gcp * psi(jl,i) / pdp(jl,i)
187     .             - hi(jl,i) * ci(jl,i+1)  )  / denom
188 
189      ai(jl,i) = -bi(jl,i) / denom
190c     -------------------
191        enddo
192      enddo
193
194
195c-------------------------------------------------------------------------
196c        2.5   
197c              -------------
198c
199
200c couche 1
201c -------
202      do jl = 1 , kdlon
203        newpcolc(jl,1) = ci(jl,1)
204      enddo
205
206
207      do i = 2 , nlaylte
208        do jl = 1 , kdlon
209           newpcolc(jl,i) = ci(jl,i) + ai(jl,i) * newpcolc(jl,i-1)
210        enddo
211      enddo
212
213
214
215c-------------------------------------------------------------------------
216      RETURN
217      END
Note: See TracBrowser for help on using the repository browser.