/[lmdze]/trunk/Sources/dyn3d/vlyqs.f
ViewVC logotype

Annotation of /trunk/Sources/dyn3d/vlyqs.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 157 - (hide annotations)
Mon Jul 20 16:01:49 2015 UTC (8 years, 10 months ago) by guez
File size: 6595 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 guez 43 SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat)
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     ! qsat est un argument de sortie pour le s-pg ....
10     !
11     !
12     ! --------------------------------------------------------------------
13     !
14 guez 139 use comconst
15 guez 43 use dimens_m
16 guez 66 use disvert_m
17 guez 57 use conf_gcm_m
18 guez 139 use comgeom, only: aire, apoln, apols
19 guez 43 USE nr_util, ONLY : pi
20 guez 139 USE dynetat0_m, only: rlonv, rlonu
21     use paramet_m
22    
23 guez 43 IMPLICIT NONE
24     !
25     !
26     ! Arguments:
27     ! ----------
28     REAL masse(ip1jmp1,llm),pente_max
29     REAL masse_adv_v( ip1jm,llm)
30     REAL q(ip1jmp1,llm)
31     REAL qsat(ip1jmp1,llm)
32     !
33     ! Local
34     ! ---------
35     !
36     INTEGER i,ij,l
37     !
38     REAL airej2,airejjm,airescb(iim),airesch(iim)
39     REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
40     REAL adyqv(ip1jm),dyqmax(ip1jmp1)
41     REAL qbyv(ip1jm,llm)
42    
43     REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
44 guez 157 Logical first
45     SAVE first
46 guez 43
47     REAL convpn,convps,convmpn,convmps
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 guez 157 DATA first/.true./
57 guez 43
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    
77    
78     DO l = 1, llm
79     !
80     ! --------------------------------
81     ! CALCUL EN LATITUDE
82     ! --------------------------------
83    
84     ! On commence par calculer la valeur du traceur moyenne sur le premier cercle
85     ! de latitude autour du pole (qpns pour le pole nord et qpsn pour
86     ! le pole nord) qui sera utilisee pour evaluer les pentes au pole.
87    
88     DO i = 1, iim
89     airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
90     airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
91     ENDDO
92     qpns = SSUM( iim, airescb ,1 ) / airej2
93     qpsn = SSUM( iim, airesch ,1 ) / airejjm
94    
95     ! calcul des pentes aux points v
96    
97     DO ij=1,ip1jm
98     dyqv(ij)=q(ij,l)-q(ij+iip1,l)
99     adyqv(ij)=abs(dyqv(ij))
100     ENDDO
101    
102     ! calcul des pentes aux points scalaires
103    
104     DO ij=iip2,ip1jm
105     dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
106     dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
107     dyqmax(ij)=pente_max*dyqmax(ij)
108     ENDDO
109    
110     ! calcul des pentes aux poles
111    
112     DO ij=1,iip1
113     dyq(ij,l)=qpns-q(ij+iip1,l)
114     dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
115     ENDDO
116    
117     ! filtrage de la derivee
118     dyn1=0.
119     dys1=0.
120     dyn2=0.
121     dys2=0.
122     DO ij=1,iim
123     dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
124     dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
125     dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
126     dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
127     ENDDO
128     DO ij=1,iip1
129     dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
130     dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
131     ENDDO
132    
133     ! calcul des pentes limites aux poles
134    
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    
150     ! calcul des pentes limitees
151    
152     DO ij=iip2,ip1jm
153     IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
154     dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
155     ELSE
156     dyq(ij,l)=0.
157     ENDIF
158     ENDDO
159    
160     ENDDO
161    
162     DO l=1,llm
163     DO ij=1,ip1jm
164     IF( masse_adv_v(ij,l).GT.0. ) THEN
165     qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l ) + &
166     dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
167     ELSE
168     qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) * &
169     0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
170     ENDIF
171     qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
172     ENDDO
173     ENDDO
174    
175    
176     DO l=1,llm
177     DO ij=iip2,ip1jm
178     newmasse=masse(ij,l) &
179     +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
180     q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l)) &
181     /newmasse
182     masse(ij,l)=newmasse
183     ENDDO
184     !.-. ancienne version
185     convpn=SSUM(iim,qbyv(1,l),1)/apoln
186     convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
187     DO ij = 1,iip1
188     newmasse=masse(ij,l)+convmpn*aire(ij)
189     q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/ &
190     newmasse
191     masse(ij,l)=newmasse
192     ENDDO
193     convps = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
194     convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
195     DO ij = ip1jm+1,ip1jmp1
196     newmasse=masse(ij,l)+convmps*aire(ij)
197     q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/ &
198     newmasse
199     masse(ij,l)=newmasse
200     ENDDO
201     !.-. fin ancienne version
202    
203     !._. nouvelle version
204     ! convpn=SSUM(iim,qbyv(1,l),1)
205     ! convmpn=ssum(iim,masse_adv_v(1,l),1)
206     ! oldmasse=ssum(iim,masse(1,l),1)
207     ! newmasse=oldmasse+convmpn
208     ! newq=(q(1,l)*oldmasse+convpn)/newmasse
209     ! newmasse=newmasse/apoln
210     ! DO ij = 1,iip1
211     ! q(ij,l)=newq
212     ! masse(ij,l)=newmasse*aire(ij)
213     ! ENDDO
214     ! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
215     ! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
216     ! oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
217     ! newmasse=oldmasse+convmps
218     ! newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
219     ! newmasse=newmasse/apols
220     ! DO ij = ip1jm+1,ip1jmp1
221     ! q(ij,l)=newq
222     ! masse(ij,l)=newmasse*aire(ij)
223     ! ENDDO
224     !._. fin nouvelle version
225     ENDDO
226    
227     RETURN
228     END

  ViewVC Help
Powered by ViewVC 1.1.21