My Project
Functions/Subroutines
adv_s.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine adv_s
 

Function/Subroutine Documentation

◆ adv_s()

subroutine adv_s ( )

Definition at line 46 of file adv_s.f90.

46 
47  !------------------------------------------------------------------------------|
48 
49  USE all_vars
50  USE mod_utils
51  USE bcs
52  USE mod_obcs
53  USE mod_par
54  USE mod_wd
55  USE mod_spherical
56  USE mod_northpole
57 
58 
59  IMPLICIT NONE
60  REAL(SP), DIMENSION(0:MT,KB) :: XFLUX,XFLUX_ADV
61  REAL(SP), DIMENSION(0:MT) :: PUPX,PUPY,PVPX,PVPY
62  REAL(SP), DIMENSION(0:MT) :: PSPX,PSPY,PSPXD,PSPYD,VISCOFF
63  REAL(SP), DIMENSION(3*(NT),KBM1) :: DTIJ
64  REAL(SP), DIMENSION(3*(NT),KBM1) :: UVN
65  REAL(SP) :: UTMP,VTMP,SITAI,FFD,FF1 !,X11,Y11,X22,Y22,X33,Y33 !,TMP1,TMP2 !,XI,YI
66  REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2,UN
67  REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF,EXFLUX,TEMP,STPOINT
68  REAL(SP) :: FACT,FM1
69  INTEGER :: I,I1,I2,IA,IB,J,J1,J2,K,JTMP,JJ,II
70  REAL(SP) :: S1MIN, S1MAX, S2MIN, S2MAX
71 
72 !!$# if defined (SPHERICAL)
73 !!$ REAL(DP) :: TY,TXPI,TYPI
74 !!$ REAL(DP) :: XTMP1,XTMP
75 !!$ REAL(DP) :: X1_DP,Y1_DP,X2_DP,Y2_DP,XII,YII
76 !!$ REAL(DP) :: X11_TMP,Y11_TMP,X33_TMP,Y33_TMP
77 !!$ REAL(DP) :: VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP
78 !!$ REAL(DP) :: TXPI_TMP,TYPI_TMP
79 !!$# endif
80 
81 
82 
83 
84 
85  IF(dbg_set(dbg_sbr)) WRITE(ipt,*)"Start: adv_s"
86  !------------------------------------------------------------------------------!
87 
88 
89  SELECT CASE(horizontal_mixing_type)
90  CASE ('closure')
91  fact = 1.0_sp
92  fm1 = 0.0_sp
93  CASE('constant')
94  fact = 0.0_sp
95  fm1 = 1.0_sp
96  CASE DEFAULT
97  CALL fatal_error("UNKNOW HORIZONTAL MIXING TYPE:",&
98  & trim(horizontal_mixing_type) )
99  END SELECT
100 
101  !
102  !--Initialize Fluxes-----------------------------------------------------------!
103  !
104  xflux = 0.0_sp
105  xflux_adv = 0.0_sp
106 
107  !
108  !--Loop Over Control Volume Sub-Edges And Calculate Normal Velocity------------!
109  !
110 !!# if !defined (1)
111  DO i=1,ncv
112  i1=ntrg(i)
113  ! DTIJ(I)=DT1(I1)
114  DO k=1,kbm1
115  dtij(i,k) = dt1(i1)*dz1(i1,k)
116  ! USE U,V
117  uvn(i,k) = v(i1,k)*dltxe(i) - u(i1,k)*dltye(i)
118  END DO
119  END DO
120 !!# else
121 !! DO I=1,NCV
122 !! I1=NTRG(I)
123 !! ! DTIJ(I)=DT1(I1)
124 !! DO K=1,KBM1
125 !! DTIJ(I,K) = DT1(I1)*DZ1(I1,K)
126 !! ! USE US,VS
127 !! UVN(I,K) = VS(I1,K)*DLTXE(I) - US(I1,K)*DLTYE(I)
128 !!# if defined (SEMI_IMPLICIT)
129 !! DTIJ1(I,K) = D1(I1)*DZ1(I1,K)
130 !! UVN1(I,K) = VF(I1,K)*DLTXE(I) - UF(I1,K)*DLTYE(I)
131 !!# endif
132 !! END DO
133 !! END DO
134 !!# endif
135 
136  !
137  !--Calculate the Advection and Horizontal Diffusion Terms----------------------!
138  !
139 
140  DO k=1,kbm1
141  pspx = 0.0_sp
142  pspy = 0.0_sp
143  pspxd = 0.0_sp
144  pspyd = 0.0_sp
145  DO i=1,m
146  DO j=1,ntsn(i)-1
147  i1=nbsn(i,j)
148  i2=nbsn(i,j+1)
149 
150 !J. Ge for tracer advection
151  IF(backward_advection==.false.)THEN
152 
153  IF(iswetn(i1) == 0 .AND. iswetn(i2) == 1)THEN
154  ffd=0.5_sp*(s1(i,k)+s1(i2,k)-smean1(i,k)-smean1(i2,k))
155  ff1=0.5_sp*(s1(i,k)+s1(i2,k))
156  ELSE IF(iswetn(i1) == 1 .AND. iswetn(i2) == 0)THEN
157  ffd=0.5_sp*(s1(i1,k)+s1(i,k)-smean1(i1,k)-smean1(i,k))
158  ff1=0.5_sp*(s1(i1,k)+s1(i,k))
159  ELSE IF(iswetn(i1) == 0 .AND. iswetn(i2) == 0)THEN
160  ffd=0.5_sp*(s1(i,k)+s1(i,k)-smean1(i,k)-smean1(i,k))
161  ff1=0.5_sp*(s1(i,k)+s1(i,k))
162  ELSE
163  ffd=0.5_sp*(s1(i1,k)+s1(i2,k)-smean1(i1,k)-smean1(i2,k))
164  ff1=0.5_sp*(s1(i1,k)+s1(i2,k))
165  END IF
166 
167  ELSE
168 
169  IF(backward_step==1)THEN
170  IF(iswetn(i1) == 0 .AND. iswetn(i2) == 1)THEN
171  ffd=0.5_sp*((s0(i,k)+s1(i,k))*0.5+(s0(i2,k)+s1(i2,k))*0.5-smean1(i,k)-smean1(i2,k))
172  ff1=0.5_sp*((s0(i,k)+s1(i,k))*0.5+(s0(i2,k)+s1(i2,k))*0.5)
173  ELSE IF(iswetn(i1) == 1 .AND. iswetn(i2) == 0)THEN
174  ffd=0.5_sp*((s0(i1,k)+s1(i1,k))*0.5+(s0(i,k)+s1(i,k))*0.5-smean1(i1,k)-smean1(i,k))
175  ff1=0.5_sp*((s0(i1,k)+s1(i1,k))*0.5+(s0(i,k)+s1(i,k))*0.5)
176  ELSE IF(iswetn(i1) == 0 .AND. iswetn(i2) == 0)THEN
177  ffd=0.5_sp*((s0(i,k)+s1(i,k))*0.5+(s0(i,k)+s1(i,k))*0.5-smean1(i,k)-smean1(i,k))
178  ff1=0.5_sp*((s0(i,k)+s1(i,k))*0.5+(s0(i,k)+s1(i,k))*0.5)
179  ELSE
180  ffd=0.5_sp*((s0(i1,k)+s1(i1,k))*0.5+(s0(i2,k)+s1(i2,k))*0.5-smean1(i1,k)-smean1(i2,k))
181  ff1=0.5_sp*((s0(i1,k)+s1(i1,k))*0.5+(s0(i2,k)+s1(i2,k))*0.5)
182  END IF
183  ELSEIF(backward_step==2)THEN
184  IF(iswetn(i1) == 0 .AND. iswetn(i2) == 1)THEN
185  ffd=0.5_sp*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp+(s2(i2,k)+s0(i2,k)+s1(i2,k))/3.0_sp-smean1(i,k)-smean1(i2,k))
186  ff1=0.5_sp*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp+(s2(i2,k)+s0(i2,k)+s1(i2,k))/3.0_sp)
187  ELSE IF(iswetn(i1) == 1 .AND. iswetn(i2) == 0)THEN
188  ffd=0.5_sp*((s2(i1,k)+s0(i1,k)+s1(i1,k))/3.0_sp+(s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp-smean1(i1,k)-smean1(i,k))
189  ff1=0.5_sp*((s2(i1,k)+s0(i1,k)+s1(i1,k))/3.0_sp+(s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp)
190  ELSE IF(iswetn(i1) == 0 .AND. iswetn(i2) == 0)THEN
191  ffd=0.5_sp*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp+(s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp-smean1(i,k)-smean1(i,k))
192  ff1=0.5_sp*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp+(s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp)
193  ELSE
194  ffd=0.5_sp*((s2(i1,k)+s0(i1,k)+s1(i1,k))/3.0_sp+(s2(i2,k)+s0(i2,k)+s1(i2,k))/3.0_sp-smean1(i1,k)-smean1(i2,k))
195  ff1=0.5_sp*((s2(i1,k)+s0(i1,k)+s1(i1,k))/3.0_sp+(s2(i2,k)+s0(i2,k)+s1(i2,k))/3.0_sp)
196  END IF
197  ENDIF
198 
199  END IF
200 !J. Ge for tracer advection
201 
202 
203 !!$# if defined (SPHERICAL)
204 !!$ XTMP = VX(I2)*TPI-VX(I1)*TPI
205 !!$ XTMP1 = VX(I2)-VX(I1)
206 !!$ IF(XTMP1 > 180.0_SP)THEN
207 !!$ XTMP = -360.0_SP*TPI+XTMP
208 !!$ ELSE IF(XTMP1 < -180.0_SP)THEN
209 !!$ XTMP = 360.0_SP*TPI+XTMP
210 !!$ END IF
211 !!$ TXPI=XTMP*COS(DEG2RAD*VY(I))
212 !!$ TYPI=(VY(I1)-VY(I2))*TPI
213 !!$ ! ERROR HERE
214 !!$!# if defined (NORTHPOLE)
215 !!$ IF(NODE_NORTHAREA(I) == 1)THEN
216 !!$ VX1_TMP = REARTH * COS(VY(I1)*DEG2RAD) * COS(VX(I1)*DEG2RAD) &
217 !!$ * 2._SP /(1._SP+sin(VY(I1)*DEG2RAD))
218 !!$ VY1_TMP = REARTH * COS(VY(I1)*DEG2RAD) * SIN(VX(I1)*DEG2RAD) &
219 !!$ * 2._SP /(1._SP+sin(VY(I1)*DEG2RAD))
220 !!$
221 !!$ VX2_TMP = REARTH * COS(VY(I2)*DEG2RAD) * COS(VX(I2)*DEG2RAD) &
222 !!$ * 2._SP /(1._SP+sin(VY(I2)*DEG2RAD))
223 !!$ VY2_TMP = REARTH * COS(VY(I2)*DEG2RAD) * SIN(VX(I2)*DEG2RAD) &
224 !!$ * 2._SP /(1._SP+sin(VY(I2)*DEG2RAD))
225 !!$
226 !!$ TXPI = (VX2_TMP-VX1_TMP)/(2._SP /(1._SP+sin(VY(I)*DEG2RAD)))
227 !!$ TYPI = (VY1_TMP-VY2_TMP)/(2._SP /(1._SP+sin(VY(I)*DEG2RAD)))
228 !!$ IF(I /= NODE_NORTHPOLE)THEN
229 !!$ TXPI_TMP = TYPI*COS(VX(I)*DEG2RAD)-TXPI*SIN(VX(I)*DEG2RAD)
230 !!$ TYPI_TMP = TXPI*COS(VX(I)*DEG2RAD)+TYPI*SIN(VX(I)*DEG2RAD)
231 !!$ TYPI_TMP = -TYPI_TMP
232 !!$
233 !!$ TXPI = TXPI_TMP
234 !!$ TYPI = TYPI_TMP
235 !!$ END IF
236 !!$ END IF
237 !!$# endif
238 !!$ ! END ERROR
239 !!$# else
240 !!$ PSPX(I)=PSPX(I)+FF1*(VY(I1)-VY(I2))
241 !!$ PSPY(I)=PSPY(I)+FF1*(VX(I2)-VX(I1))
242 !!$ PSPXD(I)=PSPXD(I)+FFD*(VY(I1)-VY(I2))
243 !!$ PSPYD(I)=PSPYD(I)+FFD*(VX(I2)-VX(I1))
244 !!$# endif
245 
246 
247  pspx(i)=pspx(i)+ff1*dltytrie(i,j)
248  pspy(i)=pspy(i)+ff1*dltxtrie(i,j)
249  pspxd(i)=pspxd(i)+ffd*dltytrie(i,j)
250  pspyd(i)=pspyd(i)+ffd*dltxtrie(i,j)
251 
252 
253  END DO
254 ! gather all neighboring control volumes connecting at dam node
255 
256  pspx(i)=pspx(i)/art2(i)
257  pspy(i)=pspy(i)/art2(i)
258  pspxd(i)=pspxd(i)/art2(i)
259  pspyd(i)=pspyd(i)/art2(i)
260 
261  END DO
262 
263  IF(k == kbm1)THEN
264  DO i=1,m
265  pfpxb(i) = pspx(i)
266  pfpyb(i) = pspy(i)
267  END DO
268  END IF
269 
270  DO i = 1,m
271  viscoff(i)=viscofh(i,k)
272  END DO
273 
274  IF(k == kbm1) THEN
275  ah_bottom(1:m) = (fact*viscoff(1:m) + fm1)*nn_hvc(1:m)
276  END IF
277 
278 
279  DO i=1,ncv_i
280  ia=niec(i,1)
281  ib=niec(i,2)
282 
283 
284 !!$ XI=0.5_SP*(XIJE(I,1)+XIJE(I,2))
285 !!$ YI=0.5_SP*(YIJE(I,1)+YIJE(I,2))
286 !!$# if defined (SPHERICAL)
287 !!$ X1_DP=XIJE(I,1)
288 !!$ Y1_DP=YIJE(I,1)
289 !!$ X2_DP=XIJE(I,2)
290 !!$ Y2_DP=YIJE(I,2)
291 !!$ CALL ARCC(X2_DP,Y2_DP,X1_DP,Y1_DP,XII,YII)
292 !!$ XI=XII
293 !!$ XTMP = XI*TPI-VX(IA)*TPI
294 !!$ XTMP1 = XI-VX(IA)
295 !!$ IF(XTMP1 > 180.0_SP)THEN
296 !!$ XTMP = -360.0_SP*TPI+XTMP
297 !!$ ELSE IF(XTMP1 < -180.0_SP)THEN
298 !!$ XTMP = 360.0_SP*TPI+XTMP
299 !!$ END IF
300 !!$
301 !!$ DXA=XTMP*COS(DEG2RAD*VY(IA))
302 !!$ DYA=(YI-VY(IA))*TPI
303 !!$ XTMP = XI*TPI-VX(IB)*TPI
304 !!$ XTMP1 = XI-VX(IB)
305 !!$ IF(XTMP1 > 180.0_SP)THEN
306 !!$ XTMP = -360.0_SP*TPI+XTMP
307 !!$ ELSE IF(XTMP1 < -180.0_SP)THEN
308 !!$ XTMP = 360.0_SP*TPI+XTMP
309 !!$ END IF
310 !!$
311 !!$ DXB=XTMP*COS(DEG2RAD*VY(IB))
312 !!$ DYB=(YI-VY(IB))*TPI
313 !!$# else
314 !!$ DXA=XI-VX(IA)
315 !!$ DYA=YI-VY(IA)
316 !!$ DXB=XI-VX(IB)
317 !!$ DYB=YI-VY(IB)
318 !!$# endif
319 
320 !J. Ge for tracer advection
321  IF(backward_advection==.false.)THEN
322  fij1=s1(ia,k)+dltxncve(i,1)*pspx(ia)+dltyncve(i,1)*pspy(ia)
323  fij2=s1(ib,k)+dltxncve(i,2)*pspx(ib)+dltyncve(i,2)*pspy(ib)
324  ELSE
325  IF(backward_step==1)THEN
326  fij1=(s0(ia,k)+s1(ia,k))*0.5+dltxncve(i,1)*pspx(ia)+dltyncve(i,1)*pspy(ia)
327  fij2=(s0(ib,k)+s1(ib,k))*0.5+dltxncve(i,2)*pspx(ib)+dltyncve(i,2)*pspy(ib)
328  ELSEIF(backward_step==2)THEN
329  fij1=(s2(ia,k)+s0(ia,k)+s1(ia,k))/3.0_sp+dltxncve(i,1)*pspx(ia)+dltyncve(i,1)*pspy(ia)
330  fij2=(s2(ib,k)+s0(ib,k)+s1(ib,k))/3.0_sp+dltxncve(i,2)*pspx(ib)+dltyncve(i,2)*pspy(ib)
331  ENDIF
332  ENDIF
333 !J. Ge for tracer advection
334  s1min=minval(s1(nbsn(ia,1:ntsn(ia)-1),k))
335  s1min=min(s1min, s1(ia,k))
336  s1max=maxval(s1(nbsn(ia,1:ntsn(ia)-1),k))
337  s1max=max(s1max, s1(ia,k))
338  s2min=minval(s1(nbsn(ib,1:ntsn(ib)-1),k))
339  s2min=min(s2min, s1(ib,k))
340  s2max=maxval(s1(nbsn(ib,1:ntsn(ib)-1),k))
341  s2max=max(s2max, s1(ib,k))
342  IF(fij1 < s1min) fij1=s1min
343  IF(fij1 > s1max) fij1=s1max
344  IF(fij2 < s2min) fij2=s2min
345  IF(fij2 > s2max) fij2=s2max
346 
347  un=uvn(i,k)
348 
349  !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
350 
351  ! David moved HPRNU and added HVC
352  viscof=(fact*0.5_sp*(viscoff(ia)*nn_hvc(ia)+viscoff(ib)*nn_hvc(ib)) + fm1*0.5_sp*(nn_hvc(ia)+nn_hvc(ib)))
353 
354  txx=0.5_sp*(pspxd(ia)+pspxd(ib))*viscof
355  tyy=0.5_sp*(pspyd(ia)+pspyd(ib))*viscof
356 
357  fxx=-dtij(i,k)*txx*dltye(i)
358  fyy= dtij(i,k)*tyy*dltxe(i)
359 
360  exflux=-un*dtij(i,k)* &
361  ((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp+fxx+fyy
362 
363  xflux(ia,k)=xflux(ia,k)+exflux
364  xflux(ib,k)=xflux(ib,k)-exflux
365 
366  xflux_adv(ia,k)=xflux_adv(ia,k)+(exflux-fxx-fyy)
367  xflux_adv(ib,k)=xflux_adv(ib,k)-(exflux-fxx-fyy)
368 
369 
370 
371  END DO
372 
373 
374 
375  END DO !!SIGMA LOOP
376 
377  !
378  !-Accumulate Fluxes at Boundary Nodes
379  !
380 
381  DO k=1,kbm1
382  IF(iobcn > 0) THEN
383  DO i=1,iobcn
384  i1=i_obc_n(i)
385  xflux_obc(i,k)=xflux_adv(i1,k)
386  END DO
387  END IF
388  END DO
389 
390  !--Set Boundary Conditions-For Fresh Water Flux--------------------------------!
391  !
392 
393 
394 
395 
396  !--------------------------------------------------------------------
397  ! The central difference scheme in vertical advection
398  !--------------------------------------------------------------------
399  DO i=1,m
400  IF(iswetn(i)*iswetnt(i) == 1) THEN
401 
402 
403  DO k=1, kbm1
404 
405 
406 
407 !J. Ge for tracer advection
408  IF(backward_advection==.false.)THEN
409  IF(k == 1) THEN
410  temp=-wts(i,k+1)*(s1(i,k)*dz(i,k+1)+s1(i,k+1)*dz(i,k))/ &
411  (dz(i,k)+dz(i,k+1))
412  ELSE IF(k == kbm1) THEN
413  temp= wts(i,k)*(s1(i,k)*dz(i,k-1)+s1(i,k-1)*dz(i,k))/(dz(i,k)+dz(i,k-1))
414  ELSE
415  temp= wts(i,k)*(s1(i,k)*dz(i,k-1)+s1(i,k-1)*dz(i,k))/(dz(i,k)+dz(i,k-1))-&
416  wts(i,k+1)*(s1(i,k)*dz(i,k+1)+s1(i,k+1)*dz(i,k))/(dz(i,k)+dz(i,k+1))
417  END IF
418  ELSE
419  IF(backward_step==1)THEN
420  IF(k == 1) THEN
421  temp=-wts(i,k+1)*((s0(i,k)+s1(i,k))*0.5*dz(i,k+1)+(s0(i,k+1)+s1(i,k+1))*0.5*dz(i,k))/ &
422  (dz(i,k)+dz(i,k+1))
423  ELSE IF(k == kbm1) THEN
424  temp= wts(i,k)*((s0(i,k)+s1(i,k))*0.5*dz(i,k-1)+(s0(i,k-1)+s1(i,k-1))*0.5*dz(i,k))/(dz(i,k)+dz(i,k-1))
425  ELSE
426  temp= wts(i,k)*((s0(i,k)+s1(i,k))*0.5*dz(i,k-1)+(s0(i,k-1)+s1(i,k-1))*0.5*dz(i,k))/(dz(i,k)+dz(i,k-1))-&
427  wts(i,k+1)*((s0(i,k)+s1(i,k))*0.5*dz(i,k+1)+(s0(i,k+1)+s1(i,k+1))*0.5*dz(i,k))/(dz(i,k)+dz(i,k+1))
428  END IF
429  ELSEIF(backward_step==2)THEN
430  IF(k == 1) THEN
431  temp=-wts(i,k+1)*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp*dz(i,k+1)+(s2(i,k+1)+s0(i,k+1)+s1(i,k+1))/3.0_sp*dz(i,k))/ &
432  (dz(i,k)+dz(i,k+1))
433  ELSE IF(k == kbm1) THEN
434  temp= wts(i,k)*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp*dz(i,k-1)+(s2(i,k-1)+s0(i,k-1)+s1(i,k-1))/3.0_sp*dz(i,k))/(dz(i,k)+dz(i,k-1))
435  ELSE
436  temp= wts(i,k)*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp*dz(i,k-1)+(s2(i,k-1)+s0(i,k-1)+s1(i,k-1))/3.0_sp*dz(i,k))/(dz(i,k)+dz(i,k-1))-&
437  wts(i,k+1)*((s2(i,k)+s0(i,k)+s1(i,k))/3.0_sp*dz(i,k+1)+(s2(i,k+1)+s0(i,k+1)+s1(i,k+1))/3.0_sp*dz(i,k))/(dz(i,k)+dz(i,k+1))
438  END IF
439  ENDIF
440  ENDIF
441 !J. Ge for tracer advection
442 
443 
444  IF(isonb(i) == 2) THEN
445  ! XFLUX(I,K)=TEMP*ART1(I)/DZ(I,K)
446  xflux(i,k)=temp*art1(i)
447  ELSE
448  ! XFLUX(I,K)=XFLUX(I,K)+TEMP*ART1(I)/DZ(I,K)
449  xflux(i,k)=xflux(i,k)+temp*art1(i)
450 
451  END IF
452  ENDDO
453  END IF
454  END DO !! SIGMA LOOP
455 
456  !
457  !--Set Boundary Conditions-For Fresh Water Flux--------------------------------!
458  !
459  IF(river_ts_setting == 'calculated') THEN
460  IF(river_inflow_location == 'node') THEN
461  IF(numqbc > 0) THEN
462  DO j=1,numqbc
463  jj=inodeq(j)
464  stpoint=sdis(j)
465  DO k=1,kbm1
466  ! XFLUX(JJ,K)=XFLUX(JJ,K) - QDIS(J)*VQDIST(J,K)*STPOINT/DZ(JJ,K)
467  xflux(jj,k)=xflux(jj,k) - qdis(j)*vqdist(j,k)*stpoint
468  END DO
469  END DO
470  END IF
471  ELSE IF(river_inflow_location == 'edge') THEN
472  IF(numqbc > 0) THEN
473  DO j=1,numqbc
474  j1=n_icellq(j,1)
475  j2=n_icellq(j,2)
476  stpoint=sdis(j) !!ASK LIU SHOULD THIS BE STPOINT1(J1)/STPOINT2(J2)
477  DO k=1,kbm1
478  ! XFLUX(J1,K)=XFLUX(J1,K)- &
479  ! QDIS(J)*RDISQ(J,1)*VQDIST(J,K)*STPOINT/DZ1(J1,K)
480  ! XFLUX(J2,K)=XFLUX(J2,K)- &
481  ! QDIS(J)*RDISQ(J,2)*VQDIST(J,K)*STPOINT/DZ1(J2,K)
482  xflux(j1,k)=xflux(j1,k)-qdis(j)*rdisq(j,1)*vqdist(j,k)*stpoint
483  xflux(j2,k)=xflux(j2,k)-qdis(j)*rdisq(j,2)*vqdist(j,k)*stpoint
484  END DO
485  END DO
486  END IF
487  END IF
488  END IF
489 !---------------------------------------------------------------------
490 
491 
492  ! APPLY GROUND WATER SALINITY FORCING
493  IF(groundwater_on .and. groundwater_salt_on)THEN
494  DO i=1,m
495  xflux(i,kbm1)=xflux(i,kbm1)-bfwdis(i)*bfwslt(i)
496  END DO
497  ELSEIF(groundwater_on) THEN
498  DO i=1,m
499 ! XFLUX(I,KBM1)=XFLUX(I,KBM1)-BFWDIS(I)*S1(I,KBM1)/DZ(I,KBM1)
500 !J. Ge for tracer advection
501  IF(backward_advection==.false.)THEN
502  xflux(i,kbm1)=xflux(i,kbm1)-bfwdis(i)*s1(i,kbm1)
503  ELSE
504  IF(backward_step==1)THEN
505  xflux(i,kbm1)=xflux(i,kbm1)-bfwdis(i)*(s0(i,kbm1)+s1(i,kbm1))*0.5
506  ELSEIF(backward_step==2)THEN
507  xflux(i,kbm1)=xflux(i,kbm1)-bfwdis(i)*(s2(i,kbm1)+s0(i,kbm1)+s1(i,kbm1))/3.0_sp
508  ENDIF
509  ENDIF
510 !J. Ge for tracer advection
511  END DO
512  END IF
513 
514 
515  !--Update Salinity-------------------------------------------------------------!
516  !
517  DO i=1,m
518  IF(iswetn(i)*iswetnt(i) == 1 )THEN
519  DO k=1,kbm1
520  ! SF1(I,K)=(S1(I,K)-XFLUX(I,K)/ART1(I)*(DTI/DT(I)))*(DT(I)/DTFA(I))
521 
522  sf1(i,k)=(s1(i,k)-xflux(i,k)/art1(i)*(dti/(dt(i)*dz(i,k))))*(dt(i)/dtfa(i))
523 
524  END DO
525 
526  ELSE ! IF NOT WET S STAYS THE SAME
527  DO k=1,kbm1
528  sf1(i,k)=s1(i,k)
529  END DO
530  END IF
531  END DO
532 
533 
534  IF(dbg_set(dbg_sbr)) WRITE(ipt,*)"End: adv_s"
535 
all_vars::ntsn
integer, dimension(:), allocatable, target ntsn
Definition: mod_main.f90:1023
mod_utils
Definition: mod_utils.f90:40
all_vars::viscofh
real(sp), dimension(:,:), allocatable, target viscofh
Definition: mod_main.f90:1359
all_vars::dtfa
real(sp), dimension(:), allocatable, target dtfa
Definition: mod_main.f90:1124
mod_spherical
Definition: mod_spherical.f90:45
all_vars::v
real(sp), dimension(:,:), allocatable, target v
Definition: mod_main.f90:1269
all_vars::vqdist
real(sp), dimension(:,:), allocatable, target vqdist
Definition: mod_main.f90:1217
mod_utils::dbg_set
logical function dbg_set(vrb)
Definition: mod_utils.f90:182
all_vars::pfpxb
real(sp), dimension(:), allocatable, target pfpxb
Definition: mod_main.f90:1342
all_vars::dltxncve
real(sp), dimension(:,:), allocatable, target dltxncve
Definition: mod_main.f90:1059
all_vars::art1
real(sp), dimension(:), allocatable, target art1
Definition: mod_main.f90:1010
all_vars::qdis
real(sp), dimension(:), allocatable, target qdis
Definition: mod_main.f90:1220
mod_obcs::xflux_obc
real(sp), dimension(:,:), allocatable xflux_obc
Definition: mod_obcs.f90:113
all_vars::dltytrie
real(sp), dimension(:,:), allocatable, target dltytrie
Definition: mod_main.f90:1063
all_vars::dltyncve
real(sp), dimension(:,:), allocatable, target dltyncve
Definition: mod_main.f90:1060
mod_par
Definition: mod_par.f90:40
all_vars::pfpyb
real(sp), dimension(:), allocatable, target pfpyb
Definition: mod_main.f90:1343
all_vars::art2
real(sp), dimension(:), allocatable, target art2
Definition: mod_main.f90:1011
all_vars::ntrg
integer, dimension(:), allocatable, target ntrg
Definition: mod_main.f90:1033
all_vars::u
real(sp), dimension(:,:), allocatable, target u
Definition: mod_main.f90:1268
all_vars::s1
real(sp), dimension(:,:), allocatable, target s1
Definition: mod_main.f90:1308
mod_northpole
Definition: mod_northpole.f90:40
all_vars::sdis
real(sp), dimension(:), allocatable, target sdis
Definition: mod_main.f90:1225
all_vars::niec
integer, dimension(:,:), allocatable, target niec
Definition: mod_main.f90:1032
all_vars::smean1
real(sp), dimension(:,:), allocatable, target smean1
Definition: mod_main.f90:1319
all_vars::rdisq
real(sp), dimension(:,:), allocatable, target rdisq
Definition: mod_main.f90:1227
all_vars::dltye
real(sp), dimension(:), allocatable, target dltye
Definition: mod_main.f90:1051
all_vars::bfwslt
real(sp), dimension(:), allocatable, target bfwslt
Definition: mod_main.f90:1201
bcs::iobcn
integer iobcn
Definition: mod_main.f90:1777
all_vars::nn_hvc
real(sp), dimension(:), allocatable nn_hvc
Definition: mod_main.f90:1303
bcs::i_obc_n
integer, dimension(:), allocatable i_obc_n
Definition: mod_main.f90:1779
all_vars::sf1
real(sp), dimension(:,:), allocatable, target sf1
Definition: mod_main.f90:1311
all_vars::bfwdis
real(sp), dimension(:), allocatable, target bfwdis
Definition: mod_main.f90:1196
all_vars::n_icellq
integer, dimension(:,:), allocatable, target n_icellq
Definition: mod_main.f90:1216
all_vars::dz
real(sp), dimension(:,:), allocatable, target dz
Definition: mod_main.f90:1092
all_vars
Definition: mod_main.f90:959
mod_wd::iswetnt
integer, dimension(:), allocatable iswetnt
Definition: mod_wd.f90:53
all_vars::dltxtrie
real(sp), dimension(:,:), allocatable, target dltxtrie
Definition: mod_main.f90:1064
all_vars::dt1
real(sp), dimension(:), allocatable, target dt1
Definition: mod_main.f90:1117
mod_wd
Definition: mod_wd.f90:40
mod_obcs
Definition: mod_obcs.f90:62
mod_utils::fatal_error
subroutine fatal_error(ER1, ER2, ER3, ER4)
Definition: mod_utils.f90:230
all_vars::dz1
real(sp), dimension(:,:), allocatable, target dz1
Definition: mod_main.f90:1096
all_vars::nbsn
integer, dimension(:,:), allocatable, target nbsn
Definition: mod_main.f90:1030
all_vars::wts
real(sp), dimension(:,:), allocatable, target wts
Definition: mod_main.f90:1321
all_vars::ah_bottom
real(sp), dimension(:), allocatable, target ah_bottom
Definition: mod_main.f90:1345
all_vars::s2
real(sp), dimension(:,:), allocatable, target s2
Definition: mod_main.f90:1316
all_vars::dltxe
real(sp), dimension(:), allocatable, target dltxe
Definition: mod_main.f90:1050
mod_utils::dbg_sbr
integer, parameter dbg_sbr
Definition: mod_utils.f90:69
all_vars::s0
real(sp), dimension(:,:), allocatable, target s0
Definition: mod_main.f90:1315
all_vars::inodeq
integer, dimension(:), allocatable, target inodeq
Definition: mod_main.f90:1214
all_vars::isonb
integer, dimension(:), allocatable, target isonb
Definition: mod_main.f90:1024
all_vars::dt
real(sp), dimension(:), allocatable, target dt
Definition: mod_main.f90:1133
mod_wd::iswetn
integer, dimension(:), allocatable iswetn
Definition: mod_wd.f90:51
Here is the call graph for this function:
Here is the caller graph for this function:
Generated by   doxygen 1.8.15