Context Navigation

lwi.F @ 1112

Last change on this file since 1112 was 1047, checked in by emillour, 11 years ago

Mars GCM:

IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
Also updated EOF routines, everything is now in eofdump_mod.F90
Removed unused routine lectfux.F (in dyn3d)

File size: 5.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format