/[lmdze]/trunk/dyn3d/Vlsplt/vly.f
ViewVC logotype

Contents of /trunk/dyn3d/Vlsplt/vly.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 157 - (show annotations)
Mon Jul 20 16:01:49 2015 UTC (8 years, 10 months ago) by guez
Original Path: trunk/Sources/dyn3d/Vlsplt/vly.f
File size: 6963 byte(s)
Just encapsulated SUBROUTINE vlsplt in a module and cleaned it.

In procedure vlx, local variables dxqu and adxqu only need indices
iip2:ip1jm. Otherwise, just cleaned vlx.

Procedures dynredem0 and dynredem1 no longer have argument fichnom,
they just operate on a file named "restart.nc". The programming
guideline here is that gcm should not be more complex than it needs by
itself, other programs (ce0l etc.) just have to adapt to gcm. So ce0l
now creates files "restart.nc" and "restartphy.nc".

In order to facilitate decentralizing the writing of "restartphy.nc",
created a procedure phyredem0 out of phyredem. phyredem0 creates the
NetCDF header of "restartphy.nc" while phyredem writes the NetCDF
variables. As the global attribute itau_phy needs to be filled in
phyredem0, at the beginnig of the run, we must compute its value
instead of just using itap. So we have a dummy argument lmt_pas of
phyredem0. Also, the ncid of "startphy.nc" is upgraded from local
variable of phyetat0 to dummy argument. phyetat0 no longer closes
"startphy.nc".

Following the same decentralizing objective, the ncid of "restart.nc"
is upgraded from local variable of dynredem0 to module variable of
dynredem0_m. "restart.nc" is not closed at the end of dynredem0 nor
opened at the beginning of dynredem1.

In procedure etat0, instead of creating many vectors of size klon
which will be filled with zeroes, just create one array null_array.

In procedure phytrac, instead of writing trs(: 1) to a text file,
write it to "restartphy.nc" (following LMDZ). This is better because
now trs(: 1) is next to its coordinates. We can write to
"restartphy.nc" from phytrac directly, and not add trs(: 1) to the
long list of variables in physiq, thanks to the decentralizing of
"restartphy.nc".

In procedure phyetat0, we no longer write to standard output the
minimum and maximum values of read arrays. It is ok to check input and
abort on invalid values but just printing statistics on input seems too
much useless computation and out of place clutter.

1 SUBROUTINE vly(q,pente_max,masse,masse_adv_v)
2 !
3 ! Auteurs: P.Le Van, F.Hourdin, F.Forget
4 !
5 ! ***********************************************************
6 ! Shema d'advection " pseudo amont " .
7 ! *************************************************************
8 ! q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
9 ! dq sont des arguments de sortie pour le s-pg ....
10
11 use comconst
12 use comgeom, only: aire
13 use conf_gcm_m
14 use dimens_m
15 use disvert_m
16 USE dynetat0_m, only: rlonv, rlonu
17 USE nr_util, ONLY : pi
18 use paramet_m
19
20 IMPLICIT NONE
21 !
22 !
23 !
24 ! Arguments:
25 ! ----------
26 REAL masse(ip1jmp1,llm)
27 real, intent(in):: pente_max
28 REAL masse_adv_v( ip1jm,llm)
29 REAL q(ip1jmp1,llm)
30 !
31 ! Local
32 ! ---------
33 !
34 INTEGER i,ij,l
35 !
36 REAL airej2,airejjm,airescb(iim),airesch(iim)
37 REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
38 REAL adyqv(ip1jm),dyqmax(ip1jmp1)
39 REAL qbyv(ip1jm,llm)
40
41 REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
42 ! REAL newq,oldmasse
43 Logical first
44 SAVE first
45
46 REAL convpn,convps,convmpn,convmps
47 real massepn,masseps,qpn,qps
48 REAL sinlon(iip1),sinlondlon(iip1)
49 REAL coslon(iip1),coslondlon(iip1)
50 SAVE sinlon,coslon,sinlondlon,coslondlon
51 SAVE airej2,airejjm
52 !
53 !
54 REAL SSUM
55
56 DATA first/.true./
57
58 IF(first) THEN
59 PRINT*,'Shema Amont nouveau appele dans Vanleer '
60 first=.false.
61 do i=2,iip1
62 coslon(i)=cos(rlonv(i))
63 sinlon(i)=sin(rlonv(i))
64 coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
65 sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
66 ENDDO
67 coslon(1)=coslon(iip1)
68 coslondlon(1)=coslondlon(iip1)
69 sinlon(1)=sinlon(iip1)
70 sinlondlon(1)=sinlondlon(iip1)
71 airej2 = SSUM( iim, aire(iip2), 1 )
72 airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
73 ENDIF
74
75 !
76 DO l = 1, llm
77 !
78 ! --------------------------------
79 ! CALCUL EN LATITUDE
80 ! --------------------------------
81
82 ! On commence par calculer la valeur du traceur moyenne sur le
83 ! premier cercle
84 ! de latitude autour du pole (qpns pour le pole nord et qpsn pour
85 ! le pole nord) qui sera utilisee pour evaluer les pentes au pole.
86
87 DO i = 1, iim
88 airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
89 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
90 ENDDO
91 qpns = SSUM( iim, airescb ,1 ) / airej2
92 qpsn = SSUM( iim, airesch ,1 ) / airejjm
93
94 ! calcul des pentes aux points v
95
96 DO ij=1,ip1jm
97 dyqv(ij)=q(ij,l)-q(ij+iip1,l)
98 adyqv(ij)=abs(dyqv(ij))
99 ENDDO
100
101 ! calcul des pentes aux points scalaires
102
103 DO ij=iip2,ip1jm
104 dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
105 dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
106 dyqmax(ij)=pente_max*dyqmax(ij)
107 ENDDO
108
109 ! calcul des pentes aux poles
110
111 DO ij=1,iip1
112 dyq(ij,l)=qpns-q(ij+iip1,l)
113 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
114 ENDDO
115
116 ! filtrage de la derivee
117 dyn1=0.
118 dys1=0.
119 dyn2=0.
120 dys2=0.
121 DO ij=1,iim
122 dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
123 dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
124 dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
125 dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
126 ENDDO
127 DO ij=1,iip1
128 dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
129 dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
130 ENDDO
131
132 ! calcul des pentes limites aux poles
133
134 goto 8888
135 fn=1.
136 fs=1.
137 DO ij=1,iim
138 IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
139 fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
140 ENDIF
141 IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
142 fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
143 ENDIF
144 ENDDO
145 DO ij=1,iip1
146 dyq(ij,l)=fn*dyq(ij,l)
147 dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
148 ENDDO
149 8888 continue
150 DO ij=1,iip1
151 dyq(ij,l)=0.
152 dyq(ip1jm+ij,l)=0.
153 ENDDO
154
155 ! calcul des pentes limitees
156
157 DO ij=iip2,ip1jm
158 IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
159 dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
160 ELSE
161 dyq(ij,l)=0.
162 ENDIF
163 ENDDO
164
165 ENDDO
166
167 DO l=1,llm
168 DO ij=1,ip1jm
169 IF(masse_adv_v(ij,l).gt.0) THEN
170 qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)* &
171 & 0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
172 ELSE
173 qbyv(ij,l)=q(ij,l)-dyq(ij,l)* &
174 & 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
175 ENDIF
176 qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
177 ENDDO
178 ENDDO
179
180
181 DO l=1,llm
182 DO ij=iip2,ip1jm
183 newmasse=masse(ij,l) &
184 & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
185 q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l)) &
186 & /newmasse
187 masse(ij,l)=newmasse
188 ENDDO
189 !.-. ancienne version
190 ! convpn=SSUM(iim,qbyv(1,l),1)/apoln
191 ! convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
192
193 convpn=SSUM(iim,qbyv(1,l),1)
194 convmpn=ssum(iim,masse_adv_v(1,l),1)
195 massepn=ssum(iim,masse(1,l),1)
196 qpn=0.
197 do ij=1,iim
198 qpn=qpn+masse(ij,l)*q(ij,l)
199 enddo
200 qpn=(qpn+convpn)/(massepn+convmpn)
201 do ij=1,iip1
202 q(ij,l)=qpn
203 enddo
204
205 ! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
206 ! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
207
208 convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
209 convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
210 masseps=ssum(iim, masse(ip1jm+1,l),1)
211 qps=0.
212 do ij = ip1jm+1,ip1jmp1-1
213 qps=qps+masse(ij,l)*q(ij,l)
214 enddo
215 qps=(qps+convps)/(masseps+convmps)
216 do ij=ip1jm+1,ip1jmp1
217 q(ij,l)=qps
218 enddo
219
220 !.-. fin ancienne version
221
222 !._. nouvelle version
223 ! convpn=SSUM(iim,qbyv(1,l),1)
224 ! convmpn=ssum(iim,masse_adv_v(1,l),1)
225 ! oldmasse=ssum(iim,masse(1,l),1)
226 ! newmasse=oldmasse+convmpn
227 ! newq=(q(1,l)*oldmasse+convpn)/newmasse
228 ! newmasse=newmasse/apoln
229 ! DO ij = 1,iip1
230 ! q(ij,l)=newq
231 ! masse(ij,l)=newmasse*aire(ij)
232 ! ENDDO
233 ! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
234 ! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
235 ! oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
236 ! newmasse=oldmasse+convmps
237 ! newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
238 ! newmasse=newmasse/apols
239 ! DO ij = ip1jm+1,ip1jmp1
240 ! q(ij,l)=newq
241 ! masse(ij,l)=newmasse*aire(ij)
242 ! ENDDO
243 !._. fin nouvelle version
244 ENDDO
245
246 RETURN
247 END

  ViewVC Help
Powered by ViewVC 1.1.21