My Project
mod_northpole.f90
Go to the documentation of this file.
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 !/===========================================================================/
13 ! Copyright (c) 2007, The University of Massachusetts Dartmouth
14 ! Produced at the School of Marine Science & Technology
15 ! Marine Ecosystem Dynamics Modeling group
16 ! All rights reserved.
17 !
18 ! FVCOM has been developed by the joint UMASSD-WHOI research team. For
19 ! details of authorship and attribution of credit please see the FVCOM
20 ! technical manual or contact the MEDM group.
21 !
22 !
23 ! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu
24 ! The full copyright notice is contained in the file COPYRIGHT located in the
25 ! root directory of the FVCOM code. This original header must be maintained
26 ! in all distributed versions.
27 !
28 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
29 ! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
30 ! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
31 ! PURPOSE ARE DISCLAIMED.
32 !
33 !/---------------------------------------------------------------------------/
34 ! CVS VERSION INFORMATION
35 ! $Id$
36 ! $Name$
37 ! $Revision$
38 !/===========================================================================/
39 
40 MODULE mod_northpole
41  USE all_vars
42  USE mod_spherical
43  USE mod_par
44  USE mod_utils
45 
46  IMPLICIT NONE
47  SAVE
48  INTEGER :: node_northpole !Node index at the north pole point
49  INTEGER :: mp,np,npe,npcv
50  INTEGER, ALLOCATABLE :: node_northarea(:)
51  INTEGER, ALLOCATABLE :: cell_northarea(:)
52  INTEGER, ALLOCATABLE :: npedge_lst(:)
53  INTEGER, ALLOCATABLE :: ncedge_lst(:)
54  INTEGER, ALLOCATABLE :: mp_lst(:),np_lst(:)
55  REAL(dp), ALLOCATABLE :: a1u_xy(:,:),a2u_xy(:,:)
56  REAL(dp), ALLOCATABLE :: aw0_xy(:,:),awx_xy(:,:),awy_xy(:,:)
57 
58 
59  CONTAINS
60 !==============================================================================|
61  SUBROUTINE find_northpole
62  IMPLICIT NONE
63  INTEGER :: I,ITMP,NODE_NORTHPOLE_GL,IERR,NDOM
64  INTEGER,ALLOCATABLE :: TMP(:)
65 
66  node_northpole = 0
67 ! NODE_NORTHPOLE_GL = 0
68  np = 0
69  mp = 0
70  npe = 0
71  npcv = 0
72 
73  ndom=0
74  DO i = 1,mt
75 ! IF(ABS(VY(I)-90.0_SP) .LE. (NEAREST(90.0_SP,1.0_sp) - 90.0_SP))THEN
76  IF(abs(vy(i)-90.0_sp) .LE. 10e-4)THEN
77  node_northpole = i
78  ndom=ndom+1
79  END IF
80  END DO
81 
82  IF (ndom .GT. 1) CALL fatal_error('Found more than one north pole node in the grid',&
83  &'This should never happen, TGE should crash first?')
84 
85 
86  ! SINCE SPHERICAL NOW ALWAYS LOOKS FOR A NORTH POLE, DO NOT CALL
87  ! FATAL_ERROR IF WE DON'T FIND A NORTH POLE!
88 
89  node_northpole_gl = ngid_x(node_northpole)
90 
91 
92 
93  if(dbg_set(dbg_log)) THEN
94  WRITE(ipt,*) "! //////////////////////////////////////////////"
95  IF(ndom .eq.0) THEN
96  WRITE(ipt,*) "! NO NORTH POLE FOUND; PROCEED WITH SPHERICAL MODEL!"
97  ELSE
98  WRITE(ipt,*) "! FOUND THE GLOBAL NORTH POLE NODE::", node_northpole_gl
99  IF(ndom .gt. 1)THEN
100  WRITE(ipt,*) "THE NORTH POLE IS ON A PROCESSOR BOUNDARY"
101  WRITE(ipt,*) "Be afraid, be very afraid...."
102  CALL fatal_error("THE NORTH POLE IS ON A PROCESSOR BOUNDARY")
103  END IF
104  END IF
105  WRITE(ipt,*) "! //////////////////////////////////////////////"
106  end if
107 
108 ! ! ALL OTHER PROCESSORS ESCAPE HERE
109 ! IF (NODE_NORTHPOLE .EQ. 0) RETURN
110 
111  ALLOCATE(node_northarea(0:mt)); node_northarea = 0
112  ALLOCATE(cell_northarea(0:nt)); cell_northarea = 0
113 
114  ! ALL OTHER PROCESSORS ESCAPE HERE
115  IF (node_northpole .EQ. 0) RETURN
116 
117 
118  itmp = node_northpole
119  ALLOCATE(tmp(mt)); tmp = 0
120  mp = 0
121  DO i=1,ntsn(itmp)-1
122  mp = mp + 1
123  tmp(mp) = nbsn(itmp,i)
124  node_northarea(nbsn(itmp,i)) = 1
125  END DO
126  mp = mp + 1
127  tmp(mp) = itmp
128  node_northarea(itmp) = 1
129 
130  ALLOCATE(mp_lst(mp))
131  mp_lst(1:mp) = tmp(1:mp)
132  DEALLOCATE(tmp)
133 
134  ALLOCATE(tmp(nt)); tmp = 0
135  np = 0
136  DO i=1,ntve(itmp)
137  np = np + 1
138  tmp(np) = nbve(itmp,i)
139  cell_northarea(nbve(itmp,i)) = 1
140  END DO
141 
142  ALLOCATE(np_lst(np))
143  np_lst(1:np) = tmp(1:np)
144  DEALLOCATE(tmp)
145 
146 
147  RETURN
148  END SUBROUTINE find_northpole
149 !==============================================================================|
150 
151 !==============================================================================|
152 
153  SUBROUTINE find_cellside
154 
155  IMPLICIT NONE
156  INTEGER :: I,IA,IB
157  INTEGER, ALLOCATABLE :: TEMP(:)
158 
159  ! ALL OTHER PROCESSORS ESCAPE HERE
160  IF (node_northpole .EQ. 0) RETURN
161 
162 
163  ALLOCATE(temp(ne)); temp = zero
164  npe = 0
165 
166  DO i=1,ne
167  ia = iec(i,1)
168  ib = iec(i,2)
169  IF(cell_northarea(ia) == 1 .OR. cell_northarea(ib) == 1)THEN
170  npe = npe + 1
171  temp(npe) = i
172  END IF
173  END DO
174 
175  ALLOCATE(npedge_lst(npe))
176  npedge_lst(1:npe) = temp(1:npe)
177  DEALLOCATE(temp)
178 
179  ALLOCATE(temp(ncv)); temp = zero
180  npcv = 0
181 
182  DO i=1,ncv
183  ia = niec(i,1)
184  ib = niec(i,2)
185  IF(ia == node_northpole .OR. ib == node_northpole)THEN
186  npcv = npcv + 1
187  temp(npcv) = i
188  END IF
189  END DO
190 
191  ALLOCATE(ncedge_lst(npcv))
192  ncedge_lst(1:npcv) = temp(1:npcv)
193  DEALLOCATE(temp)
194 
195  RETURN
196  END SUBROUTINE find_cellside
197 
198 !==============================================================================|
199 
200 !==============================================================================|
201  SUBROUTINE advave_edge_xy(XFLUX,YFLUX,IFCETA)
202 
203  USE mod_obcs
204  IMPLICIT NONE
205  INTEGER :: I,J,K,IA,IB,J1,J2,K1,K2,K3,I1,I2,II
206  REAL(SP) :: DIJ,ELIJ,XIJ,YIJ,UIJ,VIJ
207  REAL(SP) :: COFA1,COFA2,COFA3,COFA4,COFA5,COFA6,COFA7,COFA8
208  REAL(SP) :: XADV,YADV,TXXIJ,TYYIJ,TXYIJ
209  REAL(SP) :: VISCOF,VISCOF1,VISCOF2,TEMP
210  REAL(SP) :: XFLUX(0:NT),YFLUX(0:NT)
211  REAL(SP) :: FACT,FM1,ISWETTMP
212 
213  REAL(SP) :: TPA,TPB
214 
215  REAL(SP) :: UIJ1_TMP,VIJ1_TMP,UIJ2_TMP,VIJ2_TMP,TXXIJ_TMP,TYYIJ_TMP
216  REAL(SP) :: XADV_TMP,YADV_TMP,PSTX_TMP,PSTY_TMP
217  REAL(SP) :: UIJ_TMP,VIJ_TMP,UN_TMP
218  REAL(SP) :: DLTXC_TMP,DLTYC_TMP
219  REAL(SP) :: VX1_TMP,VX2_TMP,VY1_TMP,VY2_TMP
220  REAL(SP) :: UAIA,VAIA,UAIB,VAIB,UAK1,VAK1,UAK2,VAK2,UAK3,VAK3
221  REAL(SP) :: XIJC_TMP,YIJC_TMP,XCIA_TMP,YCIA_TMP,XCIB_TMP,YCIB_TMP
222  REAL(SP) :: XIJ_TMP,YIJ_TMP
223 
224  REAL(SP) :: IFCETA
225 
226  REAL(SP) :: UIJ1,VIJ1,UIJ2,VIJ2,FXX,FYY
227 !------------------------------------------------------------------------------!
228 
229  ! ALL OTHER PROCESSORS ESCAPE HERE
230  IF (node_northpole .EQ. 0) RETURN
231 
232  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: advave_edge_XY"
233 
234  SELECT CASE(horizontal_mixing_type)
235  CASE ('closure')
236  fact = 1.0_sp
237  fm1 = 0.0_sp
238  CASE('constant')
239  fact = 0.0_sp
240  fm1 = 1.0_sp
241  CASE DEFAULT
242  CALL fatal_error("UNKNOW HORIZONTAL MIXING TYPE:",&
243  & trim(horizontal_mixing_type) )
244  END SELECT
245 
246 !
247 !--Initialize Variables--------------------------------------------------------|
248 !
249  DO ii=1,npe
250  i=npedge_lst(ii)
251  ia=iec(i,1)
252  ib=iec(i,2)
253  IF(cell_northarea(ia) == 1)THEN
254  xflux(ia) = 0.0_sp
255  yflux(ia) = 0.0_sp
256  pstx(ia) = 0.0_sp
257  psty(ia) = 0.0_sp
258  END IF
259  IF(cell_northarea(ib) == 1)THEN
260  xflux(ib) = 0.0_sp
261  yflux(ib) = 0.0_sp
262  pstx(ib) = 0.0_sp
263  psty(ib) = 0.0_sp
264  END IF
265  END DO
266 !
267 !-------------------------ACCUMULATE FLUX OVER ELEMENT EDGES-------------------!
268 !
269 
270  DO ii=1,npe
271  i=npedge_lst(ii)
272  ia=iec(i,1)
273  ib=iec(i,2)
274  j1=ienode(i,1)
275  j2=ienode(i,2)
276 
277  dij=0.5_sp*(d(j1)+d(j2))
278  elij=0.5_sp*(el(j1)+el(j2))
279 
280 ! ggao 0103-2007 consider the sea surface pressure change
281 ! change end
282 
283 
284 
285 ! FLUX FROM LEFT
286  k1=nbe(ia,1)
287  k2=nbe(ia,2)
288  k3=nbe(ia,3)
289 
290  uaia = -va(ia)*cos(xc(ia)*deg2rad)-ua(ia)*sin(xc(ia)*deg2rad)
291  vaia = -va(ia)*sin(xc(ia)*deg2rad)+ua(ia)*cos(xc(ia)*deg2rad)
292  uak1 = -va(k1)*cos(xc(k1)*deg2rad)-ua(k1)*sin(xc(k1)*deg2rad)
293  vak1 = -va(k1)*sin(xc(k1)*deg2rad)+ua(k1)*cos(xc(k1)*deg2rad)
294  uak2 = -va(k2)*cos(xc(k2)*deg2rad)-ua(k2)*sin(xc(k2)*deg2rad)
295  vak2 = -va(k2)*sin(xc(k2)*deg2rad)+ua(k2)*cos(xc(k2)*deg2rad)
296  uak3 = -va(k3)*cos(xc(k3)*deg2rad)-ua(k3)*sin(xc(k3)*deg2rad)
297  vak3 = -va(k3)*sin(xc(k3)*deg2rad)+ua(k3)*cos(xc(k3)*deg2rad)
298 
299  cofa1=a1u_xy(ia,1)*uaia+a1u_xy(ia,2)*uak1 &
300  +a1u_xy(ia,3)*uak2+a1u_xy(ia,4)*uak3
301  cofa2=a2u_xy(ia,1)*uaia+a2u_xy(ia,2)*uak1 &
302  +a2u_xy(ia,3)*uak2+a2u_xy(ia,4)*uak3
303  cofa5=a1u_xy(ia,1)*vaia+a1u_xy(ia,2)*vak1 &
304  +a1u_xy(ia,3)*vak2+a1u_xy(ia,4)*vak3
305  cofa6=a2u_xy(ia,1)*vaia+a2u_xy(ia,2)*vak1 &
306  +a2u_xy(ia,3)*vak2+a2u_xy(ia,4)*vak3
307 
308  xijc_tmp = rearth * cos(yijc(i)*deg2rad) * cos(xijc(i)*deg2rad) &
309  * 2._sp /(1._sp+sin(yijc(i)*deg2rad))
310  yijc_tmp = rearth * cos(yijc(i)*deg2rad) * sin(xijc(i)*deg2rad) &
311  * 2._sp /(1._sp+sin(yijc(i)*deg2rad))
312  xcia_tmp = rearth * cos(yc(ia)*deg2rad) * cos(xc(ia)*deg2rad) &
313  * 2._sp /(1._sp+sin(yc(ia)*deg2rad))
314  ycia_tmp = rearth * cos(yc(ia)*deg2rad) * sin(xc(ia)*deg2rad) &
315  * 2._sp /(1._sp+sin(yc(ia)*deg2rad))
316 
317  xij_tmp = xijc_tmp-xcia_tmp
318  yij_tmp = yijc_tmp-ycia_tmp
319 
320  uij1=uaia+cofa1*xij_tmp+cofa2*yij_tmp
321  vij1=vaia+cofa5*xij_tmp+cofa6*xij_tmp
322 
323 ! FLUX FROM RIGHT
324  k1=nbe(ib,1)
325  k2=nbe(ib,2)
326  k3=nbe(ib,3)
327 
328  uaib = -va(ib)*cos(xc(ib)*deg2rad)-ua(ib)*sin(xc(ib)*deg2rad)
329  vaib = -va(ib)*sin(xc(ib)*deg2rad)+ua(ib)*cos(xc(ib)*deg2rad)
330  uak1 = -va(k1)*cos(xc(k1)*deg2rad)-ua(k1)*sin(xc(k1)*deg2rad)
331  vak1 = -va(k1)*sin(xc(k1)*deg2rad)+ua(k1)*cos(xc(k1)*deg2rad)
332  uak2 = -va(k2)*cos(xc(k2)*deg2rad)-ua(k2)*sin(xc(k2)*deg2rad)
333  vak2 = -va(k2)*sin(xc(k2)*deg2rad)+ua(k2)*cos(xc(k2)*deg2rad)
334  uak3 = -va(k3)*cos(xc(k3)*deg2rad)-ua(k3)*sin(xc(k3)*deg2rad)
335  vak3 = -va(k3)*sin(xc(k3)*deg2rad)+ua(k3)*cos(xc(k3)*deg2rad)
336 
337  cofa3=a1u_xy(ib,1)*uaib+a1u_xy(ib,2)*uak1 &
338  +a1u_xy(ib,3)*uak2+a1u_xy(ib,4)*uak3
339  cofa4=a2u_xy(ib,1)*uaib+a2u_xy(ib,2)*uak1 &
340  +a2u_xy(ib,3)*uak2+a2u_xy(ib,4)*uak3
341  cofa7=a1u_xy(ib,1)*vaib+a1u_xy(ib,2)*vak1 &
342  +a1u_xy(ib,3)*vak2+a1u_xy(ib,4)*vak3
343  cofa8=a2u_xy(ib,1)*vaib+a2u_xy(ib,2)*vak1 &
344  +a2u_xy(ib,3)*vak2+a2u_xy(ib,4)*vak3
345 
346  xcib_tmp = rearth * cos(yc(ib)*deg2rad) * cos(xc(ib)*deg2rad) &
347  * 2._sp /(1._sp+sin(yc(ib)*deg2rad))
348  ycib_tmp = rearth * cos(yc(ib)*deg2rad) * sin(xc(ib)*deg2rad) &
349  * 2._sp /(1._sp+sin(yc(ib)*deg2rad))
350 
351  xij_tmp = xijc_tmp-xcib_tmp
352  yij_tmp = yijc_tmp-ycib_tmp
353  uij2=uaib+cofa3*xij_tmp+cofa4*yij_tmp
354  vij2=vaib+cofa7*xij_tmp+cofa8*yij_tmp
355 
356 ! VISCOSITY COEFFICIENT
357  viscof1=art(ia)*sqrt(cofa1**2+cofa6**2+0.5_sp*(cofa2+cofa5)**2)
358  viscof2=art(ib)*sqrt(cofa3**2+cofa8**2+0.5_sp*(cofa4+cofa7)**2)
359 
360  !VISCOF=HORCON*(FACT*0.5_SP*(VISCOF1+VISCOF2)/HPRNU + FM1)
361  ! David moved HPRNU and added VHC
362  viscof=(fact*0.5_sp*(viscof1*cc_hvc(ia)+viscof2*cc_hvc(ib)) + fm1*0.5_sp*(cc_hvc(ia)+cc_hvc(ib)))/hprnu
363 
364 
365  vx1_tmp = rearth * cos(vy(ienode(i,1))*deg2rad) * cos(vx(ienode(i,1))*deg2rad) &
366  * 2._sp /(1._sp+sin(vy(ienode(i,1))*deg2rad))
367  vy1_tmp = rearth * cos(vy(ienode(i,1))*deg2rad) * sin(vx(ienode(i,1))*deg2rad) &
368  * 2._sp /(1._sp+sin(vy(ienode(i,1))*deg2rad))
369 
370  vx2_tmp = rearth * cos(vy(ienode(i,2))*deg2rad) * cos(vx(ienode(i,2))*deg2rad) &
371  * 2._sp /(1._sp+sin(vy(ienode(i,2))*deg2rad))
372  vy2_tmp = rearth * cos(vy(ienode(i,2))*deg2rad) * sin(vx(ienode(i,2))*deg2rad) &
373  * 2._sp /(1._sp+sin(vy(ienode(i,2))*deg2rad))
374 
375  dltxc_tmp = vx2_tmp-vx1_tmp
376  dltyc_tmp = vy2_tmp-vy1_tmp
377 
378 ! SHEAR STRESSES
379  txxij=(cofa1+cofa3)*viscof
380  tyyij=(cofa6+cofa8)*viscof
381  txyij=0.5_sp*(cofa2+cofa4+cofa5+cofa7)*viscof
382  fxx=dij*(txxij*dltyc_tmp-txyij*dltxc_tmp)
383  fyy=dij*(txyij*dltyc_tmp-tyyij*dltxc_tmp)
384 
385  uij1_tmp = uij1
386  vij1_tmp = vij1
387  uij2_tmp = uij2
388  vij2_tmp = vij2
389 
390  uij_tmp=0.5_sp*(uij1+uij2)
391  vij_tmp=0.5_sp*(vij1+vij2)
392  un_tmp=-uij_tmp*dltyc_tmp + vij_tmp*dltxc_tmp
393 
394 ! ADD CONVECTIVE AND VISCOUS FLUXES
395 
396 ! ACCUMULATE FLUX
397  xadv_tmp=dij*un_tmp*&
398  ((1.0_sp-sign(1.0_sp,un_tmp))*uij2_tmp &
399  +(1.0_sp+sign(1.0_sp,un_tmp))*uij1_tmp)*0.5_sp
400  yadv_tmp=dij*un_tmp* &
401  ((1.0_sp-sign(1.0_sp,un_tmp))*vij2_tmp &
402  +(1.0_sp+sign(1.0_sp,un_tmp))*vij1_tmp)*0.5_sp
403 
404  IF(cell_northarea(ia) == 1 .AND. cell_northarea(ib) == 1)THEN
405  xflux(ia)=xflux(ia)+(xadv_tmp+fxx*epor(ia))*(1.0_sp-isbc(i))*iucp(ia)
406  yflux(ia)=yflux(ia)+(yadv_tmp+fyy*epor(ia))*(1.0_sp-isbc(i))*iucp(ia)
407  xflux(ib)=xflux(ib)-(xadv_tmp+fxx*epor(ib))*(1.0_sp-isbc(i))*iucp(ib)
408  yflux(ib)=yflux(ib)-(yadv_tmp+fyy*epor(ib))*(1.0_sp-isbc(i))*iucp(ib)
409  ELSE IF(cell_northarea(ia) == 1 .AND. cell_northarea(ib) /= 1)THEN
410  xflux(ia)=xflux(ia)+(xadv_tmp+fxx*epor(ia))*(1.0_sp-isbc(i))*iucp(ia)
411  yflux(ia)=yflux(ia)+(yadv_tmp+fyy*epor(ia))*(1.0_sp-isbc(i))*iucp(ia)
412  ELSE IF(cell_northarea(ib) == 1 .AND. cell_northarea(ia) /= 1)THEN
413  xflux(ib)=xflux(ib)-(xadv_tmp+fxx*epor(ib))*(1.0_sp-isbc(i))*iucp(ib)
414  yflux(ib)=yflux(ib)-(yadv_tmp+fyy*epor(ib))*(1.0_sp-isbc(i))*iucp(ib)
415  END IF
416 
417 
418 ! ACCUMULATE BAROTROPIC FLUX
419 
420 ! VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD) &
421 ! * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))
422 ! VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD) &
423 ! * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))
424 
425 ! VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD) &
426 ! * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))
427 ! VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD) &
428 ! * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))
429 
430 ! DLTXC_TMP = VX2_TMP-VX1_TMP
431 ! DLTYC_TMP = VY2_TMP-VY1_TMP
432 
433 
434  IF(cell_northarea(ia) == 1 .AND. cell_northarea(ib) == 1)THEN
435  pstx(ia)=pstx(ia)-grav_e(ia)*d1(ia)*elij*dltyc_tmp/(2._sp /(1._sp+sin(yc(ia)*deg2rad)))
436  psty(ia)=psty(ia)+grav_e(ia)*d1(ia)*elij*dltxc_tmp/(2._sp /(1._sp+sin(yc(ia)*deg2rad)))
437  pstx(ib)=pstx(ib)+grav_e(ib)*d1(ib)*elij*dltyc_tmp/(2._sp /(1._sp+sin(yc(ib)*deg2rad)))
438  psty(ib)=psty(ib)-grav_e(ib)*d1(ib)*elij*dltxc_tmp/(2._sp /(1._sp+sin(yc(ib)*deg2rad)))
439  ELSE IF(cell_northarea(ia) == 1 .AND. cell_northarea(ib) /= 1)THEN
440  pstx(ia)=pstx(ia)-grav_e(ia)*d1(ia)*elij*dltyc_tmp/(2._sp /(1._sp+sin(yc(ia)*deg2rad)))
441  psty(ia)=psty(ia)+grav_e(ia)*d1(ia)*elij*dltxc_tmp/(2._sp /(1._sp+sin(yc(ia)*deg2rad)))
442  ELSE IF(cell_northarea(ib) == 1 .AND. cell_northarea(ia) /= 1)THEN
443  pstx(ib)=pstx(ib)+grav_e(ib)*d1(ib)*elij*dltyc_tmp/(2._sp /(1._sp+sin(yc(ib)*deg2rad)))
444  psty(ib)=psty(ib)-grav_e(ib)*d1(ib)*elij*dltxc_tmp/(2._sp /(1._sp+sin(yc(ib)*deg2rad)))
445  END IF
446 
447  END DO
448 
449  if(dbg_set(dbg_sbr)) write(ipt,*) "End: advave_edge_XY"
450 
451  RETURN
452  END SUBROUTINE advave_edge_xy
453 !==============================================================================|
454 
455 
456 
457 !==============================================================================|
458  SUBROUTINE advection_edge_xy(XFLUX,YFLUX)
459 
460  IMPLICIT NONE
461 ! REAL(SP), INTENT(OUT), DIMENSION(0:NT,KB) :: XFLUX,YFLUX
462  REAL(SP) :: XFLUX(0:NT,KB),YFLUX(0:NT,KB)
463  REAL(SP) :: DIJ
464  REAL(SP) :: COFA1,COFA2,COFA3,COFA4,COFA5,COFA6,COFA7,COFA8
465  REAL(SP) :: XADV,YADV,TXXIJ,TYYIJ,TXYIJ,UN
466  REAL(SP) :: VISCOF,VISCOF1,VISCOF2,TEMP,TPA,TPB
467  REAL(SP) :: XIJA,YIJA,XIJB,YIJB,UIJ,VIJ
468  REAL(SP) :: FACT,FM1
469  INTEGER :: I,IA,IB,J1,J2,K1,K2,K3,K4,K5,K6,K,II,J,I1,I2
470 
471  REAL(SP) :: UIJ1_TMP,VIJ1_TMP,UIJ2_TMP,VIJ2_TMP,TXXIJ_TMP,TYYIJ_TMP
472  REAL(SP) :: XADV_TMP,YADV_TMP
473  REAL(SP) :: UIJ_TMP,VIJ_TMP,UN_TMP
474  REAL(SP) :: DLTXC_TMP,DLTYC_TMP
475  REAL(SP) :: VX1_TMP,VX2_TMP,VY1_TMP,VY2_TMP
476  REAL(SP) :: UIA,VIA,UIB,VIB,UK1,VK1,UK2,VK2,UK3,VK3,UK4,VK4,UK5,VK5,UK6,VK6
477  REAL(SP) :: XIJC_TMP,YIJC_TMP,XCIA_TMP,YCIA_TMP,XCIB_TMP,YCIB_TMP
478 
479  REAL(SP) :: UIJ1,VIJ1,UIJ2,VIJ2,FXX,FYY
480 !------------------------------------------------------------------------------|
481  ! ALL OTHER PROCESSORS ESCAPE HERE
482  IF (node_northpole .EQ. 0) RETURN
483 
484  if(dbg_set(dbg_sbr)) write(ipt,*) "Start: advection_edge_XY"
485 
486  SELECT CASE(horizontal_mixing_type)
487  CASE ('closure')
488  fact = 1.0_sp
489  fm1 = 0.0_sp
490  CASE('constant')
491  fact = 0.0_sp
492  fm1 = 1.0_sp
493  CASE DEFAULT
494  CALL fatal_error("UNKNOW HORIZONTAL MIXING TYPE:",&
495  & trim(horizontal_mixing_type) )
496  END SELECT
497 
498 !
499 !--Initialize Variables--------------------------------------------------------|
500 !
501  DO k = 1,kbm1
502  DO ii=1,npe
503  i=npedge_lst(ii)
504  ia=iec(i,1)
505  ib=iec(i,2)
506  IF(cell_northarea(ia) == 1)THEN
507  xflux(ia,k) = 0.0_sp
508  yflux(ia,k) = 0.0_sp
509  END IF
510  IF(cell_northarea(ib) == 1)THEN
511  xflux(ib,k) = 0.0_sp
512  yflux(ib,k) = 0.0_sp
513  END IF
514  END DO
515  END DO
516 
517 !
518 !--Loop Over Edges and Accumulate Fluxes-For Each Element----------------------|
519 !
520 
521 
522  DO ii=1,npe
523  i=npedge_lst(ii)
524  ia=iec(i,1)
525  ib=iec(i,2)
526  j1=ienode(i,1)
527  j2=ienode(i,2)
528 ! DIJ= 0.5_SP*(DT(J1)+DT(J2))
529 
530  k1=nbe(ia,1)
531  k2=nbe(ia,2)
532  k3=nbe(ia,3)
533  k4=nbe(ib,1)
534  k5=nbe(ib,2)
535  k6=nbe(ib,3)
536 
537  xijc_tmp = rearth * cos(yijc(i)*deg2rad) * cos(xijc(i)*deg2rad) &
538  * 2._sp /(1._sp+sin(yijc(i)*deg2rad))
539  yijc_tmp = rearth * cos(yijc(i)*deg2rad) * sin(xijc(i)*deg2rad) &
540  * 2._sp /(1._sp+sin(yijc(i)*deg2rad))
541  xcia_tmp = rearth * cos(yc(ia)*deg2rad) * cos(xc(ia)*deg2rad) &
542  * 2._sp /(1._sp+sin(yc(ia)*deg2rad))
543  ycia_tmp = rearth * cos(yc(ia)*deg2rad) * sin(xc(ia)*deg2rad) &
544  * 2._sp /(1._sp+sin(yc(ia)*deg2rad))
545  xcib_tmp = rearth * cos(yc(ib)*deg2rad) * cos(xc(ib)*deg2rad) &
546  * 2._sp /(1._sp+sin(yc(ib)*deg2rad))
547  ycib_tmp = rearth * cos(yc(ib)*deg2rad) * sin(xc(ib)*deg2rad) &
548  * 2._sp /(1._sp+sin(yc(ib)*deg2rad))
549 
550  xija = xijc_tmp-xcia_tmp
551  yija = yijc_tmp-ycia_tmp
552  xijb = xijc_tmp-xcib_tmp
553  yijb = yijc_tmp-ycib_tmp
554 
555  DO k=1,kbm1
556 
557  dij= 0.5_sp*(dt(j1)*dz(j1,k)+dt(j2)*dz(j2,k))
558 
559  uia = -v(ia,k)*cos(xc(ia)*deg2rad)-u(ia,k)*sin(xc(ia)*deg2rad)
560  via = -v(ia,k)*sin(xc(ia)*deg2rad)+u(ia,k)*cos(xc(ia)*deg2rad)
561  uib = -v(ib,k)*cos(xc(ib)*deg2rad)-u(ib,k)*sin(xc(ib)*deg2rad)
562  vib = -v(ib,k)*sin(xc(ib)*deg2rad)+u(ib,k)*cos(xc(ib)*deg2rad)
563  uk1 = -v(k1,k)*cos(xc(k1)*deg2rad)-u(k1,k)*sin(xc(k1)*deg2rad)
564  vk1 = -v(k1,k)*sin(xc(k1)*deg2rad)+u(k1,k)*cos(xc(k1)*deg2rad)
565  uk2 = -v(k2,k)*cos(xc(k2)*deg2rad)-u(k2,k)*sin(xc(k2)*deg2rad)
566  vk2 = -v(k2,k)*sin(xc(k2)*deg2rad)+u(k2,k)*cos(xc(k2)*deg2rad)
567  uk3 = -v(k3,k)*cos(xc(k3)*deg2rad)-u(k3,k)*sin(xc(k3)*deg2rad)
568  vk3 = -v(k3,k)*sin(xc(k3)*deg2rad)+u(k3,k)*cos(xc(k3)*deg2rad)
569  uk4 = -v(k4,k)*cos(xc(k4)*deg2rad)-u(k4,k)*sin(xc(k4)*deg2rad)
570  vk4 = -v(k4,k)*sin(xc(k4)*deg2rad)+u(k4,k)*cos(xc(k4)*deg2rad)
571  uk5 = -v(k5,k)*cos(xc(k5)*deg2rad)-u(k5,k)*sin(xc(k5)*deg2rad)
572  vk5 = -v(k5,k)*sin(xc(k5)*deg2rad)+u(k5,k)*cos(xc(k5)*deg2rad)
573  uk6 = -v(k6,k)*cos(xc(k6)*deg2rad)-u(k6,k)*sin(xc(k6)*deg2rad)
574  vk6 = -v(k6,k)*sin(xc(k6)*deg2rad)+u(k6,k)*cos(xc(k6)*deg2rad)
575  !!FORM THE LEFT FLUX
576  cofa1=a1u_xy(ia,1)*uia+a1u_xy(ia,2)*uk1 &
577  +a1u_xy(ia,3)*uk2+a1u_xy(ia,4)*uk3
578  cofa2=a2u_xy(ia,1)*uia+a2u_xy(ia,2)*uk1 &
579  +a2u_xy(ia,3)*uk2+a2u_xy(ia,4)*uk3
580  cofa5=a1u_xy(ia,1)*via+a1u_xy(ia,2)*vk1 &
581  +a1u_xy(ia,3)*vk2+a1u_xy(ia,4)*vk3
582  cofa6=a2u_xy(ia,1)*via+a2u_xy(ia,2)*vk1 &
583  +a2u_xy(ia,3)*vk2+a2u_xy(ia,4)*vk3
584  uij1=uia+cofa1*xija+cofa2*yija
585  vij1=via+cofa5*xija+cofa6*yija
586 
587  !!FORM THE RIGHT FLUX
588  cofa3=a1u_xy(ib,1)*uib+a1u_xy(ib,2)*uk4 &
589  +a1u_xy(ib,3)*uk5+a1u_xy(ib,4)*uk6
590  cofa4=a2u_xy(ib,1)*uib+a2u_xy(ib,2)*uk4 &
591  +a2u_xy(ib,3)*uk5+a2u_xy(ib,4)*uk6
592  cofa7=a1u_xy(ib,1)*vib+a1u_xy(ib,2)*vk4 &
593  +a1u_xy(ib,3)*vk5+a1u_xy(ib,4)*vk6
594  cofa8=a2u_xy(ib,1)*vib+a2u_xy(ib,2)*vk4 &
595  +a2u_xy(ib,3)*vk5+a2u_xy(ib,4)*vk6
596  uij2=uib+cofa3*xijb+cofa4*yijb
597  vij2=vib+cofa7*xijb+cofa8*yijb
598 
599  !! VISCOSITY COEFFICIENT EDGE
600  viscof1=art(ia)*sqrt(cofa1**2+cofa6**2+0.5_sp*(cofa2+cofa5)**2)
601  viscof2=art(ib)*sqrt(cofa3**2+cofa8**2+0.5_sp*(cofa4+cofa7)**2)
602 
603 ! VISCOF = HORCON*(FACT*0.5_SP*(VISCOF1+VISCOF2)/HPRNU + FM1)
604  ! David moved HPRNU and added VHC
605  viscof=(fact*0.5_sp*(viscof1*cc_hvc(ia)+viscof2*cc_hvc(ib)) + fm1*0.5_sp*(cc_hvc(ia)+cc_hvc(ib)))/hprnu
606 
607  vx1_tmp = rearth * cos(vy(ienode(i,1))*deg2rad) * cos(vx(ienode(i,1))*deg2rad) &
608  * 2._sp /(1._sp+sin(vy(ienode(i,1))*deg2rad))
609  vy1_tmp = rearth * cos(vy(ienode(i,1))*deg2rad) * sin(vx(ienode(i,1))*deg2rad)&
610  * 2._sp /(1._sp+sin(vy(ienode(i,1))*deg2rad))
611 
612  vx2_tmp = rearth * cos(vy(ienode(i,2))*deg2rad) * cos(vx(ienode(i,2))*deg2rad)&
613  * 2._sp /(1._sp+sin(vy(ienode(i,2))*deg2rad))
614  vy2_tmp = rearth * cos(vy(ienode(i,2))*deg2rad) * sin(vx(ienode(i,2))*deg2rad)&
615  * 2._sp /(1._sp+sin(vy(ienode(i,2))*deg2rad))
616 
617  dltxc_tmp = vx2_tmp-vx1_tmp
618  dltyc_tmp = vy2_tmp-vy1_tmp
619 
620  txxij=(cofa1+cofa3)*viscof
621  tyyij=(cofa6+cofa8)*viscof
622  txyij=0.5_sp*(cofa2+cofa4+cofa5+cofa7)*viscof
623  fxx=dij*(txxij*dltyc_tmp-txyij*dltxc_tmp)
624  fyy=dij*(txyij*dltyc_tmp-tyyij*dltxc_tmp)
625 
626  uij1_tmp = uij1
627  vij1_tmp = vij1
628  uij2_tmp = uij2
629  vij2_tmp = vij2
630 
631  uij_tmp=0.5_sp*(uij1+uij2)
632  vij_tmp=0.5_sp*(vij1+vij2)
633  un_tmp=-uij_tmp*dltyc_tmp + vij_tmp*dltxc_tmp
634 
635  !!COMPUTE BOUNDARY FLUX AUGMENTERS
636  tpa = float(1-isbc(i))*epor(ia)
637  tpb = float(1-isbc(i))*epor(ib)
638 
639 
640  !!ACCUMULATE THE FLUX
641  xadv_tmp=dij*un_tmp*&
642  ((1.0_sp-sign(1.0_sp,un_tmp))*uij2_tmp &
643  +(1.0_sp+sign(1.0_sp,un_tmp))*uij1_tmp)*0.5_sp
644  yadv_tmp=dij*un_tmp* &
645  ((1.0_sp-sign(1.0_sp,un_tmp))*vij2_tmp &
646  +(1.0_sp+sign(1.0_sp,un_tmp))*vij1_tmp)*0.5_sp
647  IF(cell_northarea(ia) == 1 .AND. cell_northarea(ib) == 1)THEN
648  xflux(ia,k)=xflux(ia,k)+xadv_tmp*tpa+(fxx+3.0_sp*fxx*float(isbc(i)))*epor(ia)
649  yflux(ia,k)=yflux(ia,k)+yadv_tmp*tpa+(fyy+3.0_sp*fyy*float(isbc(i)))*epor(ia)
650  xflux(ib,k)=xflux(ib,k)-xadv_tmp*tpb-(fxx+3.0_sp*fxx*float(isbc(i)))*epor(ib)
651  yflux(ib,k)=yflux(ib,k)-yadv_tmp*tpb-(fyy+3.0_sp*fyy*float(isbc(i)))*epor(ib)
652  ELSE IF(cell_northarea(ia) == 1 .AND. cell_northarea(ib) /= 1)THEN
653  xflux(ia,k)=xflux(ia,k)+xadv_tmp*tpa+(fxx+3.0_sp*fxx*float(isbc(i)))*epor(ia)
654  yflux(ia,k)=yflux(ia,k)+yadv_tmp*tpa+(fyy+3.0_sp*fyy*float(isbc(i)))*epor(ia)
655  ELSE IF(cell_northarea(ib) == 1 .AND. cell_northarea(ia) /= 1)THEN
656  xflux(ib,k)=xflux(ib,k)-xadv_tmp*tpb-(fxx+3.0_sp*fxx*float(isbc(i)))*epor(ib)
657  yflux(ib,k)=yflux(ib,k)-yadv_tmp*tpb-(fyy+3.0_sp*fyy*float(isbc(i)))*epor(ib)
658  END IF
659  END DO
660  END DO
661 
662  if(dbg_set(dbg_sbr)) write(ipt,*) "End: advection_edge_XY"
663 
664  RETURN
665  END SUBROUTINE advection_edge_xy
666 !==============================================================================|
667 
668 !==============================================================================!
669 
670  SUBROUTINE adv_uv_edge_xy(XFLUX,YFLUX,CETA,STG,K_STG)
671 
672  IMPLICIT NONE
673  REAL(SP) :: XFLUX(0:NT,KB),YFLUX(0:NT,KB)
674  REAL(SP) :: PSTX_TM(0:NT,KB),PSTY_TM(0:NT,KB)
675  REAL(SP) :: COFA1,COFA2,COFA3,COFA4,COFA5,COFA6,COFA7,COFA8
676  REAL(SP) :: XADV,YADV,TXXIJ,TYYIJ,TXYIJ,UN
677  REAL(SP) :: VISCOF,VISCOF1,VISCOF2,TEMP,TPA,TPB
678  REAL(SP) :: XIJA,YIJA,XIJB,YIJB,UIJ,VIJ
679  REAL(SP) :: SITA,DIJ,ELIJ,TMPA,TMPB,TMP,XFLUXV,YFLUXV
680  REAL(SP) :: FACT,FM1,EXFLUX,ISWETTMP
681  INTEGER :: I,IA,IB,J1,J2,K1,K2,K3,K4,K5,K6,K,II,J,I1,I2
682 
683 ! REAL(SP) :: TY
684  REAL(SP) :: UIJ1_TMP,VIJ1_TMP,UIJ2_TMP,VIJ2_TMP,UIJ3_TMP,VIJ3_TMP
685  REAL(SP) :: U_TMP,V_TMP,UF_TMP,VF_TMP
686  REAL(SP) :: TXXIJ_TMP,TYYIJ_TMP
687  REAL(SP) :: XADV_TMP,YADV_TMP,PSTX_TMP,PSTY_TMP
688  REAL(SP) :: UIJ_TMP,VIJ_TMP,EXFLUX_TMP
689  REAL(SP) :: DLTXC_TMP,DLTYC_TMP
690  REAL(SP) :: VX1_TMP,VX2_TMP,VY1_TMP,VY2_TMP
691  REAL(SP) :: UIA,VIA,UIB,VIB,UK1,VK1,UK2,VK2,UK3,VK3,UK4,VK4,UK5,VK5,UK6,VK6
692  REAL(SP) :: XIJC_TMP,YIJC_TMP,XCIA_TMP,YCIA_TMP,XCIB_TMP,YCIB_TMP
693 
694  REAL(SP) :: CETA
695 
696  REAL(SP) :: UIJ1,VIJ1,UIJ2,VIJ2,FXX,FYY
697 
698  INTEGER :: STG, K_STG
699 !------------------------------------------------------------------------------!
700  ! ALL OTHER PROCESSORS ESCAPE HERE
701  IF (node_northpole .EQ. 0) RETURN
702 
703  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: adv_uv_edge_xy"
704 
705  SELECT CASE(horizontal_mixing_type)
706  CASE ('closure')
707  fact = 1.0_sp
708  fm1 = 0.0_sp
709  CASE('constant')
710  fact = 0.0_sp
711  fm1 = 1.0_sp
712  CASE DEFAULT
713  CALL fatal_error("UNKNOW HORIZONTAL MIXING TYPE:",&
714  & trim(horizontal_mixing_type) )
715  END SELECT
716 
717 !
718 !-----Initialize Flux Variables------------------------------------------------!
719 !
720 
721  DO k = 1,kbm1
722  DO ii=1,npe
723  i=npedge_lst(ii)
724  ia=iec(i,1)
725  ib=iec(i,2)
726  IF(cell_northarea(ia) == 1)THEN
727  xflux(ia,k) = 0.0_sp
728  yflux(ia,k) = 0.0_sp
729  pstx_tm(ia,k) = 0.0_sp
730  psty_tm(ia,k) = 0.0_sp
731  END IF
732  IF(cell_northarea(ib) == 1)THEN
733  xflux(ib,k) = 0.0_sp
734  yflux(ib,k) = 0.0_sp
735  pstx_tm(ib,k) = 0.0_sp
736  psty_tm(ib,k) = 0.0_sp
737  END IF
738  END DO
739  END DO
740 
741 !
742 !-----Loop Over Edges and Accumulate Flux--------------------------------------!
743 !
744 
745  DO ii=1,npe
746  i=npedge_lst(ii)
747  ia=iec(i,1)
748  ib=iec(i,2)
749  j1=ienode(i,1)
750  j2=ienode(i,2)
751 ! DIJ=0.5_SP*(DT(J1)+DT(J2))
752 
753 
754  elij=0.5_sp*(egf(j1)+egf(j2))
755 ! ggao 0103-2007 consider the sea surface pressure change
756 ! change end
757 
758 
759 
760  k1=nbe(ia,1)
761  k2=nbe(ia,2)
762  k3=nbe(ia,3)
763  k4=nbe(ib,1)
764  k5=nbe(ib,2)
765  k6=nbe(ib,3)
766 
767  xijc_tmp = rearth * cos(yijc(i)*deg2rad) * cos(xijc(i)*deg2rad) &
768  * 2._sp /(1._sp+sin(yijc(i)*deg2rad))
769  yijc_tmp = rearth * cos(yijc(i)*deg2rad) * sin(xijc(i)*deg2rad) &
770  * 2._sp /(1._sp+sin(yijc(i)*deg2rad))
771 
772  xcia_tmp = rearth * cos(yc(ia)*deg2rad) * cos(xc(ia)*deg2rad) &
773  * 2._sp /(1._sp+sin(yc(ia)*deg2rad))
774  ycia_tmp = rearth * cos(yc(ia)*deg2rad) * sin(xc(ia)*deg2rad) &
775  * 2._sp /(1._sp+sin(yc(ia)*deg2rad))
776  xcib_tmp = rearth * cos(yc(ib)*deg2rad) * cos(xc(ib)*deg2rad) &
777  * 2._sp /(1._sp+sin(yc(ib)*deg2rad))
778  ycib_tmp = rearth * cos(yc(ib)*deg2rad) * sin(xc(ib)*deg2rad) &
779  * 2._sp /(1._sp+sin(yc(ib)*deg2rad))
780 
781  xija = xijc_tmp-xcia_tmp
782  yija = yijc_tmp-ycia_tmp
783  xijb = xijc_tmp-xcib_tmp
784  yijb = yijc_tmp-ycib_tmp
785 
786  DO k=1,kbm1
787 
788  dij=0.5_sp*(dt(j1)*dz(j1,k)+dt(j2)*dz(j2,k))
789 
790  uia = -v(ia,k)*cos(xc(ia)*deg2rad)-u(ia,k)*sin(xc(ia)*deg2rad)
791  via = -v(ia,k)*sin(xc(ia)*deg2rad)+u(ia,k)*cos(xc(ia)*deg2rad)
792  uib = -v(ib,k)*cos(xc(ib)*deg2rad)-u(ib,k)*sin(xc(ib)*deg2rad)
793  vib = -v(ib,k)*sin(xc(ib)*deg2rad)+u(ib,k)*cos(xc(ib)*deg2rad)
794  uk1 = -v(k1,k)*cos(xc(k1)*deg2rad)-u(k1,k)*sin(xc(k1)*deg2rad)
795  vk1 = -v(k1,k)*sin(xc(k1)*deg2rad)+u(k1,k)*cos(xc(k1)*deg2rad)
796  uk2 = -v(k2,k)*cos(xc(k2)*deg2rad)-u(k2,k)*sin(xc(k2)*deg2rad)
797  vk2 = -v(k2,k)*sin(xc(k2)*deg2rad)+u(k2,k)*cos(xc(k2)*deg2rad)
798  uk3 = -v(k3,k)*cos(xc(k3)*deg2rad)-u(k3,k)*sin(xc(k3)*deg2rad)
799  vk3 = -v(k3,k)*sin(xc(k3)*deg2rad)+u(k3,k)*cos(xc(k3)*deg2rad)
800  uk4 = -v(k4,k)*cos(xc(k4)*deg2rad)-u(k4,k)*sin(xc(k4)*deg2rad)
801  vk4 = -v(k4,k)*sin(xc(k4)*deg2rad)+u(k4,k)*cos(xc(k4)*deg2rad)
802  uk5 = -v(k5,k)*cos(xc(k5)*deg2rad)-u(k5,k)*sin(xc(k5)*deg2rad)
803  vk5 = -v(k5,k)*sin(xc(k5)*deg2rad)+u(k5,k)*cos(xc(k5)*deg2rad)
804  uk6 = -v(k6,k)*cos(xc(k6)*deg2rad)-u(k6,k)*sin(xc(k6)*deg2rad)
805  vk6 = -v(k6,k)*sin(xc(k6)*deg2rad)+u(k6,k)*cos(xc(k6)*deg2rad)
806 
807  cofa1=a1u_xy(ia,1)*uia+a1u_xy(ia,2)*uk1 &
808  +a1u_xy(ia,3)*uk2+a1u_xy(ia,4)*uk3
809  cofa2=a2u_xy(ia,1)*uia+a2u_xy(ia,2)*uk1 &
810  +a2u_xy(ia,3)*uk2+a2u_xy(ia,4)*uk3
811  cofa5=a1u_xy(ia,1)*via+a1u_xy(ia,2)*vk1 &
812  +a1u_xy(ia,3)*vk2+a1u_xy(ia,4)*vk3
813  cofa6=a2u_xy(ia,1)*via+a2u_xy(ia,2)*vk1 &
814  +a2u_xy(ia,3)*vk2+a2u_xy(ia,4)*vk3
815 
816  uij1=uia+cofa1*xija+cofa2*yija
817  vij1=via+cofa5*xija+cofa6*yija
818 
819  cofa3=a1u_xy(ib,1)*uib+a1u_xy(ib,2)*uk4 &
820  +a1u_xy(ib,3)*uk5+a1u_xy(ib,4)*uk6
821  cofa4=a2u_xy(ib,1)*uib+a2u_xy(ib,2)*uk4 &
822  +a2u_xy(ib,3)*uk5+a2u_xy(ib,4)*uk6
823  cofa7=a1u_xy(ib,1)*vib+a1u_xy(ib,2)*vk4 &
824  +a1u_xy(ib,3)*vk5+a1u_xy(ib,4)*vk6
825  cofa8=a2u_xy(ib,1)*vib+a2u_xy(ib,2)*vk4 &
826  +a2u_xy(ib,3)*vk5+a2u_xy(ib,4)*vk6
827 
828  uij2=uib+cofa3*xijb+cofa4*yijb
829  vij2=vib+cofa7*xijb+cofa8*yijb
830 
831 !-------ADD THE VISCOUS TERM & ADVECTION TERM---------------------------------!
832 !
833  viscof1=art(ia)*sqrt(cofa1**2+cofa6**2+0.5_sp*(cofa2+cofa5)**2)
834  viscof2=art(ib)*sqrt(cofa3**2+cofa8**2+0.5_sp*(cofa4+cofa7)**2)
835 
836  !VISCOF=FACT*0.5_SP*HORCON*(VISCOF1+VISCOF2)/HPRNU + FM1*HORCON
837 
838  ! David moved HPRNU and added VHC
839  viscof=(fact*0.5_sp*(viscof1*cc_hvc(ia)+viscof2*cc_hvc(ib)) + fm1*0.5_sp*(cc_hvc(ia)+cc_hvc(ib)))/hprnu
840 
841  vx1_tmp = rearth * cos(vy(ienode(i,1))*deg2rad) * cos(vx(ienode(i,1))*deg2rad) &
842  * 2._sp /(1._sp+sin(vy(ienode(i,1))*deg2rad))
843  vy1_tmp = rearth * cos(vy(ienode(i,1))*deg2rad) * sin(vx(ienode(i,1))*deg2rad)&
844  * 2._sp /(1._sp+sin(vy(ienode(i,1))*deg2rad))
845 
846  vx2_tmp = rearth * cos(vy(ienode(i,2))*deg2rad) * cos(vx(ienode(i,2))*deg2rad)&
847  * 2._sp /(1._sp+sin(vy(ienode(i,2))*deg2rad))
848  vy2_tmp = rearth * cos(vy(ienode(i,2))*deg2rad) * sin(vx(ienode(i,2))*deg2rad)&
849  * 2._sp /(1._sp+sin(vy(ienode(i,2))*deg2rad))
850 
851  dltxc_tmp = vx2_tmp-vx1_tmp
852  dltyc_tmp = vy2_tmp-vy1_tmp
853 
854  txxij=(cofa1+cofa3)*viscof
855  tyyij=(cofa6+cofa8)*viscof
856  txyij=0.5_sp*(cofa2+cofa4+cofa5+cofa7)*viscof
857  fxx=dij*(txxij*dltyc_tmp-txyij*dltxc_tmp)
858  fyy=dij*(txyij*dltyc_tmp-tyyij*dltxc_tmp)
859 
860 
861  uij1_tmp = uij1
862  vij1_tmp = vij1
863  uij2_tmp = uij2
864  vij2_tmp = vij2
865 
866  uij_tmp=0.5_sp*(uij1+uij2)
867  vij_tmp=0.5_sp*(vij1+vij2)
868  exflux_tmp = dij*(-uij_tmp*dltyc_tmp + vij_tmp*dltxc_tmp)
869 
870  !!CALCULATE BOUNDARY FLUX AUGMENTERS
871  tpa = float(1-isbc(i))*epor(ia)
872  tpb = float(1-isbc(i))*epor(ib)
873 
874  !!ACCUMULATE ADVECTIVE + DIFFUSIVE + BAROTROPIC PRESSURE GRADIENT TERMS
875  xadv_tmp=exflux_tmp*&
876  ((1.0_sp-sign(1.0_sp,exflux_tmp))*uij2_tmp &
877  +(1.0_sp+sign(1.0_sp,exflux_tmp))*uij1_tmp)*0.5_sp
878  yadv_tmp=exflux_tmp* &
879  ((1.0_sp-sign(1.0_sp,exflux_tmp))*vij2_tmp &
880  +(1.0_sp+sign(1.0_sp,exflux_tmp))*vij1_tmp)*0.5_sp
881  IF(cell_northarea(ia) == 1 .AND. cell_northarea(ib) == 1)THEN
882  xflux(ia,k)=xflux(ia,k)+xadv_tmp*tpa+(fxx+3.0_sp*fxx*float(isbc(i)))*epor(ia)
883  yflux(ia,k)=yflux(ia,k)+yadv_tmp*tpa+(fyy+3.0_sp*fyy*float(isbc(i)))*epor(ia)
884  xflux(ib,k)=xflux(ib,k)-xadv_tmp*tpb-(fxx+3.0_sp*fxx*float(isbc(i)))*epor(ib)
885  yflux(ib,k)=yflux(ib,k)-yadv_tmp*tpb-(fyy+3.0_sp*fyy*float(isbc(i)))*epor(ib)
886  ELSE IF(cell_northarea(ia) == 1 .AND. cell_northarea(ib) /= 1)THEN
887  xflux(ia,k)=xflux(ia,k)+xadv_tmp*tpa+(fxx+3.0_sp*fxx*float(isbc(i)))*epor(ia)
888  yflux(ia,k)=yflux(ia,k)+yadv_tmp*tpa+(fyy+3.0_sp*fyy*float(isbc(i)))*epor(ia)
889  ELSE IF(cell_northarea(ib) == 1 .AND. cell_northarea(ia) /= 1)THEN
890  xflux(ib,k)=xflux(ib,k)-xadv_tmp*tpb-(fxx+3.0_sp*fxx*float(isbc(i)))*epor(ib)
891  yflux(ib,k)=yflux(ib,k)-yadv_tmp*tpb-(fyy+3.0_sp*fyy*float(isbc(i)))*epor(ib)
892  END IF
893 
894  vx1_tmp = rearth * cos(vy(ienode(i,1))*deg2rad) * cos(vx(ienode(i,1))*deg2rad) &
895  * 2._sp /(1._sp+sin(vy(ienode(i,1))*deg2rad))
896  vy1_tmp = rearth * cos(vy(ienode(i,1))*deg2rad) * sin(vx(ienode(i,1))*deg2rad) &
897  * 2._sp /(1._sp+sin(vy(ienode(i,1))*deg2rad))
898 
899  vx2_tmp = rearth * cos(vy(ienode(i,2))*deg2rad) * cos(vx(ienode(i,2))*deg2rad) &
900  * 2._sp /(1._sp+sin(vy(ienode(i,2))*deg2rad))
901  vy2_tmp = rearth * cos(vy(ienode(i,2))*deg2rad) * sin(vx(ienode(i,2))*deg2rad) &
902  * 2._sp /(1._sp+sin(vy(ienode(i,2))*deg2rad))
903 
904  dltxc_tmp = vx2_tmp-vx1_tmp
905  dltyc_tmp = vy2_tmp-vy1_tmp
906 
907  IF(cell_northarea(ia) == 1 .AND. cell_northarea(ib) == 1)THEN
908  pstx_tm(ia,k)=pstx_tm(ia,k)-grav_e(ia)*dt1(ia)*dz1(ia,k)*elij*dltyc_tmp/ &
909  (2._sp /(1._sp+sin(yc(ia)*deg2rad)))
910  psty_tm(ia,k)=psty_tm(ia,k)+grav_e(ia)*dt1(ia)*dz1(ia,k)*elij*dltxc_tmp/ &
911  (2._sp /(1._sp+sin(yc(ia)*deg2rad)))
912  pstx_tm(ib,k)=pstx_tm(ib,k)+grav_e(ib)*dt1(ib)*dz1(ib,k)*elij*dltyc_tmp/ &
913  (2._sp /(1._sp+sin(yc(ib)*deg2rad)))
914  psty_tm(ib,k)=psty_tm(ib,k)-grav_e(ib)*dt1(ib)*dz1(ib,k)*elij*dltxc_tmp/ &
915  (2._sp /(1._sp+sin(yc(ib)*deg2rad)))
916  ELSE IF(cell_northarea(ia) == 1 .AND. cell_northarea(ib) /= 1)THEN
917  pstx_tm(ia,k)=pstx_tm(ia,k)-grav_e(ia)*dt1(ia)*dz1(ia,k)*elij*dltyc_tmp/ &
918  (2._sp /(1._sp+sin(yc(ia)*deg2rad)))
919  psty_tm(ia,k)=psty_tm(ia,k)+grav_e(ia)*dt1(ia)*dz1(ia,k)*elij*dltxc_tmp/ &
920  (2._sp /(1._sp+sin(yc(ia)*deg2rad)))
921  ELSE IF(cell_northarea(ib) == 1 .AND. cell_northarea(ia) /= 1)THEN
922  pstx_tm(ib,k)=pstx_tm(ib,k)+grav_e(ib)*dt1(ib)*dz1(ib,k)*elij*dltyc_tmp/ &
923  (2._sp /(1._sp+sin(yc(ib)*deg2rad)))
924  psty_tm(ib,k)=psty_tm(ib,k)-grav_e(ib)*dt1(ib)*dz1(ib,k)*elij*dltxc_tmp/ &
925  (2._sp /(1._sp+sin(yc(ib)*deg2rad)))
926  END IF
927  END DO
928  END DO
929 
930 
931  DO k = 1,kbm1
932  DO ii=1,np
933  i=np_lst(ii)
934  xflux(i,k)=xflux(i,k)+pstx_tm(i,k)
935  yflux(i,k)=yflux(i,k)+psty_tm(i,k)
936  END DO
937  END DO
938 
939 !
940 !-------ADD VERTICAL CONVECTIVE FLUX, CORIOLIS TERM AND BAROCLINIC PG TERM----!
941 !
942  DO ii=1,np
943  i=np_lst(ii)
944  IF(cell_northarea(i) == 1)THEN
945  DO k=1,kbm1
946  IF(k == 1) THEN
947  uij1_tmp = -v(i,k)*cos(xc(i)*deg2rad)-u(i,k)*sin(xc(i)*deg2rad)
948  vij1_tmp = -v(i,k)*sin(xc(i)*deg2rad)+u(i,k)*cos(xc(i)*deg2rad)
949  uij2_tmp = -v(i,k+1)*cos(xc(i)*deg2rad)-u(i,k+1)*sin(xc(i)*deg2rad)
950  vij2_tmp = -v(i,k+1)*sin(xc(i)*deg2rad)+u(i,k+1)*cos(xc(i)*deg2rad)
951  xfluxv=-w(i,k+1)*(uij1_tmp*dz1(i,k+1)+uij2_tmp*dz1(i,k))/(dz1(i,k)+dz1(i,k+1))
952  yfluxv=-w(i,k+1)*(vij1_tmp*dz1(i,k+1)+vij2_tmp*dz1(i,k))/(dz1(i,k)+dz1(i,k+1))
953  ELSE IF(k == kbm1) THEN
954  uij1_tmp = -v(i,k)*cos(xc(i)*deg2rad)-u(i,k)*sin(xc(i)*deg2rad)
955  vij1_tmp = -v(i,k)*sin(xc(i)*deg2rad)+u(i,k)*cos(xc(i)*deg2rad)
956  uij2_tmp = -v(i,k-1)*cos(xc(i)*deg2rad)-u(i,k-1)*sin(xc(i)*deg2rad)
957  vij2_tmp = -v(i,k-1)*sin(xc(i)*deg2rad)+u(i,k-1)*cos(xc(i)*deg2rad)
958  xfluxv= w(i,k)*(uij1_tmp*dz1(i,k-1)+uij2_tmp*dz1(i,k))/(dz1(i,k)+dz1(i,k-1))
959  yfluxv= w(i,k)*(vij1_tmp*dz1(i,k-1)+vij2_tmp*dz1(i,k))/(dz1(i,k)+dz1(i,k-1))
960  ELSE
961  uij1_tmp = -v(i,k)*cos(xc(i)*deg2rad)-u(i,k)*sin(xc(i)*deg2rad)
962  vij1_tmp = -v(i,k)*sin(xc(i)*deg2rad)+u(i,k)*cos(xc(i)*deg2rad)
963  uij2_tmp = -v(i,k-1)*cos(xc(i)*deg2rad)-u(i,k-1)*sin(xc(i)*deg2rad)
964  vij2_tmp = -v(i,k-1)*sin(xc(i)*deg2rad)+u(i,k-1)*cos(xc(i)*deg2rad)
965  uij3_tmp = -v(i,k+1)*cos(xc(i)*deg2rad)-u(i,k+1)*sin(xc(i)*deg2rad)
966  vij3_tmp = -v(i,k+1)*sin(xc(i)*deg2rad)+u(i,k+1)*cos(xc(i)*deg2rad)
967  xfluxv= w(i,k)*(uij1_tmp*dz1(i,k-1)+uij2_tmp*dz1(i,k))/(dz1(i,k)+dz1(i,k-1))- &
968  w(i,k+1)*(uij1_tmp*dz1(i,k+1)+uij3_tmp*dz1(i,k))/(dz1(i,k)+dz1(i,k+1))
969  yfluxv= w(i,k)*(vij1_tmp*dz1(i,k-1)+vij2_tmp*dz1(i,k))/(dz1(i,k)+dz1(i,k-1))- &
970  w(i,k+1)*(vij1_tmp*dz1(i,k+1)+vij3_tmp*dz1(i,k))/(dz1(i,k)+dz1(i,k+1))
971  END IF
972  u_tmp = -v(i,k)*cos(xc(i)*deg2rad)-u(i,k)*sin(xc(i)*deg2rad)
973  v_tmp = -v(i,k)*sin(xc(i)*deg2rad)+u(i,k)*cos(xc(i)*deg2rad)
974 ! XFLUX(I,K)=XFLUX(I,K)+XFLUXV/DZ(K)*ART(I)&
975 ! +DRHOX(I,K)-COR(I)*V_TMP*DT1(I)*ART(I)
976 ! YFLUX(I,K)=YFLUX(I,K)+YFLUXV/DZ(K)*ART(I)&
977 ! +DRHOY(I,K)+COR(I)*U_TMP*DT1(I)*ART(I)
978  xflux(i,k)=xflux(i,k)+xfluxv*art(i)&
979  +drhox(i,k)-cor(i)*v_tmp*dt1(i)*dz1(i,k)*art(i)
980  yflux(i,k)=yflux(i,k)+yfluxv*art(i)&
981  +drhoy(i,k)+cor(i)*u_tmp*dt1(i)*dz1(i,k)*art(i)
982 
983  END DO
984  END IF
985  END DO
986 
987 
988  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "End: adv_uv_edge_xy"
989 
990  RETURN
991  END SUBROUTINE adv_uv_edge_xy
992 !==============================================================================!
993 
994 !==============================================================================|
995  SUBROUTINE extel_edge_xy(K,XFLUX)
996 
997  USE mod_obcs
998  IMPLICIT NONE
999  REAL(SP) :: XFLUX(0:MT)
1000  REAL(SP) :: DIJ,UIJ,VIJ
1001  INTEGER :: I,J,K,I1,IA,IB,JJ,J1,J2,II
1002 
1003  REAL(SP) :: UIJ_TMP,VIJ_TMP,EXFLUX_TMP
1004  REAL(SP) :: DLTXE_TMP,DLTYE_TMP
1005  REAL(SP) :: VX1_TMP,VX2_TMP,VY1_TMP,VY2_TMP
1006 
1007  ! ALL OTHER PROCESSORS ESCAPE HERE
1008  IF (node_northpole .EQ. 0) RETURN
1009 
1010  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: extel_edge_XY"
1011 !----------INITIALIZE FLUX ARRAY ----------------------------------------------!
1012 
1013  DO ii=1,npcv
1014  i = ncedge_lst(ii)
1015  ia = niec(i,1)
1016  ib = niec(i,2)
1017  IF(ia == node_northpole)THEN
1018  xflux(ia) = 0.0_sp
1019  END IF
1020  IF(ib == node_northpole)THEN
1021  xflux(ib) = 0.0_sp
1022  END IF
1023  END DO
1024 
1025 !---------ACCUMULATE FLUX BY LOOPING OVER CONTROL VOLUME HALF EDGES------------!
1026 
1027  DO ii=1,npcv
1028  i = ncedge_lst(ii)
1029  i1 = ntrg(i)
1030  ia = niec(i,1)
1031  ib = niec(i,2)
1032  dij = d1(i1)
1033 
1034  uij = ua(i1)
1035  vij = va(i1)
1036 
1037  IF(ia == node_northpole .OR. ib == node_northpole)THEN
1038  uij_tmp = -vij*cos(xc(i1)*deg2rad)-uij*sin(xc(i1)*deg2rad)
1039  vij_tmp = -vij*sin(xc(i1)*deg2rad)+uij*cos(xc(i1)*deg2rad)
1040 
1041  vx1_tmp = rearth * cos(yije(i,1)*deg2rad) * cos(xije(i,1)*deg2rad)&
1042  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
1043  vy1_tmp = rearth * cos(yije(i,1)*deg2rad) * sin(xije(i,1)*deg2rad)&
1044  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
1045 
1046  vx2_tmp = rearth * cos(yije(i,2)*deg2rad) * cos(xije(i,2)*deg2rad)&
1047  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
1048  vy2_tmp = rearth * cos(yije(i,2)*deg2rad) * sin(xije(i,2)*deg2rad)&
1049  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
1050 
1051  dltxe_tmp = vx2_tmp-vx1_tmp
1052  dltye_tmp = vy2_tmp-vy1_tmp
1053 
1054  exflux_tmp = dij*(-uij_tmp*dltye_tmp+vij_tmp*dltxe_tmp)
1055  END IF
1056 
1057  IF(ia == node_northpole) THEN
1058  xflux(ia) = xflux(ia)-exflux_tmp
1059  ELSE IF(ib == node_northpole)THEN
1060  xflux(ib) = xflux(ib)+exflux_tmp
1061  END IF
1062 
1063  END DO
1064 
1065  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "End: extel_edge_XY"
1066 
1067  RETURN
1068  END SUBROUTINE extel_edge_xy
1069 !==============================================================================|
1070 
1071 !==============================================================================|
1072  SUBROUTINE extuv_edge_xy(K)
1073 
1074  IMPLICIT NONE
1075  INTEGER, INTENT(IN) :: K
1076  REAL(SP), DIMENSION(0:NT) :: RESX,RESY
1077  INTEGER :: I,II
1078 
1079  REAL(SP) :: UARK_TMP,VARK_TMP,UAF_TMP,VAF_TMP,UA_TMP,VA_TMP
1080 
1081  REAL(SP) :: WUSURF2_TMP,WVSURF2_TMP,WUBOT_TMP,WVBOT_TMP
1082  ! ALL OTHER PROCESSORS ESCAPE HERE
1083  IF (node_northpole .EQ. 0) RETURN
1084 
1085  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: extuv_edge_XY"
1086 
1087 !
1088 !--ACCUMULATE RESIDUALS FOR EXTERNAL MODE EQUATIONS----------------------------|
1089 !
1090  DO ii=1,np
1091  i=np_lst(ii)
1092  ua_tmp = -va(i)*cos(xc(i)*deg2rad)-ua(i)*sin(xc(i)*deg2rad)
1093  va_tmp = -va(i)*sin(xc(i)*deg2rad)+ua(i)*cos(xc(i)*deg2rad)
1094 
1095  wusurf2_tmp = -wvsurf2(i)*cos(xc(i)*deg2rad)-wusurf2(i)*sin(xc(i)*deg2rad)
1096  wvsurf2_tmp = -wvsurf2(i)*sin(xc(i)*deg2rad)+wusurf2(i)*cos(xc(i)*deg2rad)
1097 
1098  wubot_tmp = -wvbot(i)*cos(xc(i)*deg2rad)-wubot(i)*sin(xc(i)*deg2rad)
1099  wvbot_tmp = -wvbot(i)*sin(xc(i)*deg2rad)+wubot(i)*cos(xc(i)*deg2rad)
1100 
1101  resx(i) = adx2d(i)+advua(i)+drx2d(i)+pstx(i) &
1102  -cor(i)*va_tmp*d1(i)*art(i) &
1103  -(wusurf2_tmp+wubot_tmp)*art(i)
1104  resy(i) = ady2d(i)+advva(i)+dry2d(i)+psty(i) &
1105  +cor(i)*ua_tmp*d1(i)*art(i) &
1106  -(wvsurf2_tmp+wvbot_tmp)*art(i)
1107 
1108 !
1109 !--UPDATE----------------------------------------------------------------------|
1110 !
1111  uark_tmp = -vark(i)*cos(xc(i)*deg2rad)-uark(i)*sin(xc(i)*deg2rad)
1112  vark_tmp = -vark(i)*sin(xc(i)*deg2rad)+uark(i)*cos(xc(i)*deg2rad)
1113 
1114  uaf_tmp = (uark_tmp*(h1(i)+elrk1(i)) &
1115  -alpha_rk(k)*dte*resx(i)/art(i))/(h1(i)+elf1(i))
1116  vaf_tmp = (vark_tmp*(h1(i)+elrk1(i)) &
1117  -alpha_rk(k)*dte*resy(i)/art(i))/(h1(i)+elf1(i))
1118 
1119  uaf(i) = vaf_tmp*cos(xc(i)*deg2rad)-uaf_tmp*sin(xc(i)*deg2rad)
1120  vaf(i) = uaf_tmp*cos(xc(i)*deg2rad)+vaf_tmp*sin(xc(i)*deg2rad)
1121  vaf(i) = -vaf(i)
1122 
1123  END DO
1124 
1125  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "End: extuv_edge_XY"
1126 
1127  RETURN
1128  END SUBROUTINE extuv_edge_xy
1129 !==============================================================================|
1130 
1131 !==============================================================================|
1132 
1133  SUBROUTINE vertvl_edge_xy(XFLUX,CETA)
1134 
1135 !------------------------------------------------------------------------------|
1136  IMPLICIT NONE
1137  REAL(SP) :: XFLUX(MT,KBM1)
1138  REAL(SP) :: DIJ,UIJ,VIJ,UN,EXFLUX,TMP1,DIJ1,UIJ1,VIJ1
1139  INTEGER :: I,K,IA,IB,I1 ,J,JJ,J1,J2,II
1140 
1141  REAL(SP) :: UIJ_TMP,VIJ_TMP,VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP,UIJ1_TMP,VIJ1_TMP
1142  REAL(SP) :: DLTXE_TMP,DLTYE_TMP,EXFLUX_TMP
1143 
1144  REAL(SP) :: CETA
1145 !------------------------------------------------------------------------------|
1146  ! ALL OTHER PROCESSORS ESCAPE HERE
1147  IF (node_northpole .EQ. 0) RETURN
1148 
1149  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: Vertvl_edge_XY"
1150 
1151 !----------------------INITIALIZE FLUX-----------------------------------------!
1152 
1153  DO k=1,kbm1
1154  DO ii=1,npcv
1155  i = ncedge_lst(ii)
1156  ia = niec(i,1)
1157  ib = niec(i,2)
1158  IF(ia == node_northpole)THEN
1159  xflux(ia,k) = 0.0_sp
1160  END IF
1161  IF(ib == node_northpole)THEN
1162  xflux(ib,k) = 0.0_sp
1163  END IF
1164  END DO
1165  END DO
1166 !----------------------ACCUMULATE FLUX-----------------------------------------!
1167 
1168  DO ii=1,npcv
1169  i=ncedge_lst(ii)
1170  i1=ntrg(i)
1171  ia=niec(i,1)
1172  ib=niec(i,2)
1173 ! DIJ=DT1(I1)
1174  DO k=1,kbm1
1175  dij=dt1(i1)*dz1(i1,k)
1176  uij=u(i1,k)
1177  vij=v(i1,k)
1178 
1179  IF(ia == node_northpole .OR. ib == node_northpole)THEN
1180  uij_tmp = -vij*cos(xc(i1)*deg2rad)-uij*sin(xc(i1)*deg2rad)
1181  vij_tmp = -vij*sin(xc(i1)*deg2rad)+uij*cos(xc(i1)*deg2rad)
1182 
1183  vx1_tmp = rearth * cos(yije(i,1)*deg2rad) * cos(xije(i,1)*deg2rad)&
1184  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
1185  vy1_tmp = rearth * cos(yije(i,1)*deg2rad) * sin(xije(i,1)*deg2rad)&
1186  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
1187 
1188  vx2_tmp = rearth * cos(yije(i,2)*deg2rad) * cos(xije(i,2)*deg2rad)&
1189  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
1190  vy2_tmp = rearth * cos(yije(i,2)*deg2rad) * sin(xije(i,2)*deg2rad)&
1191  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
1192 
1193  dltxe_tmp = vx2_tmp-vx1_tmp
1194  dltye_tmp = vy2_tmp-vy1_tmp
1195 
1196  exflux_tmp = dij*(-uij_tmp*dltye_tmp+vij_tmp*dltxe_tmp)
1197  END IF
1198 
1199  IF(ia == node_northpole)THEN
1200  xflux(ia,k) = xflux(ia,k)-exflux_tmp
1201  ELSE IF(ib == node_northpole)THEN
1202  xflux(ib,k) = xflux(ib,k)+exflux_tmp
1203  END IF
1204  END DO
1205  END DO
1206 
1207  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "End: Vertvl_edge_xy"
1208 
1209  RETURN
1210  END SUBROUTINE vertvl_edge_xy
1211 !==============================================================================|
1212 
1213 !==============================================================================|
1214  SUBROUTINE adv_s_xy(XFLUX,XFLUX_ADV,PSPX,PSPY,PSPXD,PSPYD,VISCOFF,K,CETA)
1215 
1216 !------------------------------------------------------------------------------|
1217 
1218  IMPLICIT NONE
1219  INTEGER, INTENT(IN) :: K
1220 
1221  REAL(SP), DIMENSION(0:MT,KB) :: XFLUX,XFLUX_ADV
1222  REAL(SP), DIMENSION(M) :: PSPX,PSPY,PSPXD,PSPYD,VISCOFF
1223 ! REAL(SP), DIMENSION(3*(NT)) :: DTIJ
1224  REAL(SP), DIMENSION(3*(NT),KBM1) :: DTIJ
1225  REAL(SP) :: XI,YI
1226  REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2
1227  REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF
1228  REAL(SP) :: FACT,FM1
1229  INTEGER :: I,I1,I2,IA,IB,J,J1,J2,JTMP,JJ,II
1230  REAL(SP) :: TXPI,TYPI
1231 
1232  REAL(SP) :: VX_TMP,VY_TMP,VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP,VX3_TMP,VY3_TMP
1233  REAL(SP) :: XI_TMP,YI_TMP,VXA_TMP,VYA_TMP,VXB_TMP,VYB_TMP
1234  REAL(SP) :: UIJ_TMP,VIJ_TMP,DLTXE_TMP,DLTYE_TMP,UVN_TMP,EXFLUX_TMP
1235  REAL(SP) :: PUPX_TMP,PUPY_TMP,PVPX_TMP,PVPY_TMP
1236  REAL(SP) :: PSPX_TMP,PSPY_TMP,PSPXD_TMP,PSPYD_TMP
1237  REAL(SP) :: U_TMP,V_TMP
1238  REAL(SP) :: X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2
1239 
1240 !! ggao edge calculation
1241  REAL(SP) :: XIJE1_TMP,YIJE1_TMP,XIJE2_TMP,YIJE2_TMP
1242  REAL(SP) :: S1MIN, S1MAX, S2MIN, S2MAX
1243 
1244  REAL(SP) :: CETA
1245 !------------------------------------------------------------------------------!
1246  ! ALL OTHER PROCESSORS ESCAPE HERE
1247  IF (node_northpole .EQ. 0) RETURN
1248 
1249  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: ADV_S_XY(K):",k
1250 
1251 
1252  SELECT CASE(horizontal_mixing_type)
1253  CASE ('closure')
1254  fact = 1.0_sp
1255  fm1 = 0.0_sp
1256  CASE('constant')
1257  fact = 0.0_sp
1258  fm1 = 1.0_sp
1259  CASE DEFAULT
1260  CALL fatal_error("UNKNOW HORIZONTAL MIXING TYPE:",&
1261  & trim(horizontal_mixing_type) )
1262  END SELECT
1263 
1264 !
1265 !--Initialize Fluxes-----------------------------------------------------------!
1266 !
1267  DO ii=1,npcv
1268  i = ncedge_lst(ii)
1269  ia = niec(i,1)
1270  ib = niec(i,2)
1271  IF(ia == node_northpole)THEN
1272  xflux(ia,k) = 0.0_sp
1273  xflux_adv(ia,k) = 0.0_sp
1274  ELSE IF(ib == node_northpole)THEN
1275  xflux(ib,k) = 0.0_sp
1276  xflux_adv(ib,k) = 0.0_sp
1277  END IF
1278  END DO
1279 
1280 !
1281 !--Loop Over Control Volume Sub-Edges And Calculate Normal Velocity------------!
1282 !
1283  DO ii=1,npcv
1284  i = ncedge_lst(ii)
1285  i1=ntrg(i)
1286  dtij(i,k)=dt1(i1)*dz1(i1,k)
1287  END DO
1288 
1289 !
1290 !--Calculate the Advection and Horizontal Diffusion Terms----------------------!
1291 !
1292  i = node_northpole
1293 
1294  IF(i==0) RETURN
1295 
1296  pupx_tmp=0.0_sp
1297  pupy_tmp=0.0_sp
1298  pvpx_tmp=0.0_sp
1299  pvpy_tmp=0.0_sp
1300 
1301  DO j=1,ntve(i)
1302  i1=nbve(i,j)
1303  jtmp=nbvt(i,j)
1304  j1=jtmp+1-(jtmp+1)/4*3
1305  j2=jtmp+2-(jtmp+2)/4*3
1306 
1307  vx_tmp = rearth * cos(vy(i)*deg2rad) * cos(vx(i)*deg2rad) &
1308  * 2._sp /(1._sp+sin(vy(i)*deg2rad))
1309  vy_tmp = rearth * cos(vy(i)*deg2rad) * sin(vx(i)*deg2rad) &
1310  * 2._sp /(1._sp+sin(vy(i)*deg2rad))
1311 
1312  vx1_tmp= rearth * cos(vy(nv(i1,j1))*deg2rad) * cos(vx(nv(i1,j1))*deg2rad) &
1313  * 2._sp /(1._sp+sin(vy(nv(i1,j1))*deg2rad))
1314  vy1_tmp= rearth * cos(vy(nv(i1,j1))*deg2rad) * sin(vx(nv(i1,j1))*deg2rad) &
1315  * 2._sp /(1._sp+sin(vy(nv(i1,j1))*deg2rad))
1316 
1317  vx2_tmp= rearth * cos(yc(i1)*deg2rad) * cos(xc(i1)*deg2rad) &
1318  * 2._sp /(1._sp+sin(yc(i1)*deg2rad))
1319  vy2_tmp= rearth * cos(yc(i1)*deg2rad) * sin(xc(i1)*deg2rad) &
1320  * 2._sp /(1._sp+sin(yc(i1)*deg2rad))
1321 
1322  vx3_tmp= rearth * cos(vy(nv(i1,j2))*deg2rad) * cos(vx(nv(i1,j2))*deg2rad) &
1323  * 2._sp /(1._sp+sin(vy(nv(i1,j2))*deg2rad))
1324  vy3_tmp= rearth * cos(vy(nv(i1,j2))*deg2rad) * sin(vx(nv(i1,j2))*deg2rad) &
1325  * 2._sp /(1._sp+sin(vy(nv(i1,j2))*deg2rad))
1326 
1327  x11=0.5_sp*(vx_tmp+vx1_tmp)
1328  y11=0.5_sp*(vy_tmp+vy1_tmp)
1329  x22=vx2_tmp
1330  y22=vx2_tmp
1331  x33=0.5_sp*(vx_tmp+vx3_tmp)
1332  y33=0.5_sp*(vy_tmp+vy3_tmp)
1333 
1334  u_tmp = -v(i1,k)*cos(xc(i1)*deg2rad)-u(i1,k)*sin(xc(i1)*deg2rad)
1335  v_tmp = -v(i1,k)*sin(xc(i1)*deg2rad)+u(i1,k)*cos(xc(i1)*deg2rad)
1336 
1337  pupx_tmp=pupx_tmp+u_tmp*(y11-y33)
1338  pupy_tmp=pupy_tmp+u_tmp*(x33-x11)
1339  pvpx_tmp=pvpx_tmp+v_tmp*(y11-y33)
1340  pvpy_tmp=pvpy_tmp+v_tmp*(x33-x11)
1341  END DO
1342 
1343  pupx_tmp=pupx_tmp/art1(i)
1344  pupy_tmp=pupy_tmp/art1(i)
1345  pvpx_tmp=pvpx_tmp/art1(i)
1346  pvpy_tmp=pvpy_tmp/art1(i)
1347  tmp1=pupx_tmp**2+pvpy_tmp**2
1348  tmp2=0.5_sp*(pupy_tmp+pvpx_tmp)**2
1349  viscoff(i)=sqrt(tmp1+tmp2)*art1(i)
1350 
1351  IF(k == kbm1) THEN
1352  ah_bottom(i) = (fact*viscoff(i) + fm1)*nn_hvc(i)
1353  END IF
1354 
1355  DO ii=1,npcv
1356  i = ncedge_lst(ii)
1357  i1=ntrg(i)
1358  ia=niec(i,1)
1359  ib=niec(i,2)
1360 
1361  IF((ia <= m .AND. ib <= m) .AND. i1 <= n)THEN
1362 ! XI=0.5_SP*(XIJE(I,1)+XIJE(I,2))
1363 ! YI=0.5_SP*(YIJE(I,1)+YIJE(I,2))
1364 !! ggao edge calculation
1365  xije1_tmp = rearth * cos(yije(i,1)*deg2rad) * cos(xije(i,1)*deg2rad) &
1366  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
1367  yije1_tmp = rearth * cos(yije(i,1)*deg2rad) * sin(xije(i,1)*deg2rad) &
1368  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
1369 
1370  xije2_tmp = rearth * cos(yije(i,2)*deg2rad) * cos(xije(i,2)*deg2rad) &
1371  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
1372  yije2_tmp = rearth * cos(yije(i,2)*deg2rad) * sin(xije(i,2)*deg2rad) &
1373  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
1374  xi_tmp =0.5_sp*(xije1_tmp+xije2_tmp)
1375  yi_tmp =0.5_sp*(yije1_tmp+yije2_tmp)
1376 
1377 
1378  IF(ia == node_northpole .OR. ib == node_northpole)THEN
1379 ! XI_TMP = REARTH * COS(YI*DEG2RAD) * COS(XI*DEG2RAD) &
1380 ! * 2._SP /(1._SP+sin(YI*DEG2RAD))
1381 ! YI_TMP = REARTH * COS(YI*DEG2RAD) * SIN(XI*DEG2RAD) &
1382 ! * 2._SP /(1._SP+sin(YI*DEG2RAD))
1383 
1384  vxa_tmp = rearth * cos(vy(ia)*deg2rad) * cos(vx(ia)*deg2rad) &
1385  * 2._sp /(1._sp+sin(vy(ia)*deg2rad))
1386  vya_tmp = rearth * cos(vy(ia)*deg2rad) * sin(vx(ia)*deg2rad) &
1387  * 2._sp /(1._sp+sin(vy(ia)*deg2rad))
1388 
1389  vxb_tmp = rearth * cos(vy(ib)*deg2rad) * cos(vx(ib)*deg2rad) &
1390  * 2._sp /(1._sp+sin(vy(ib)*deg2rad))
1391  vyb_tmp = rearth * cos(vy(ib)*deg2rad) * sin(vx(ib)*deg2rad) &
1392  * 2._sp /(1._sp+sin(vy(ib)*deg2rad))
1393 
1394 ! IF(IA == NODE_NORTHPOLE)THEN
1395  dxa=xi_tmp-vxa_tmp
1396  dya=yi_tmp-vya_tmp
1397 ! ELSE IF(IB == NODE_NORTHPOLE)THEN
1398  dxb=xi_tmp-vxb_tmp
1399  dyb=yi_tmp-vyb_tmp
1400 ! END IF
1401 ! END IF
1402 
1403  IF(ia == node_northpole)THEN
1404  pspx_tmp=-pspy(ib)*cos(vx(ib)*deg2rad)-pspx(ib)*sin(vx(ib)*deg2rad)
1405  pspy_tmp=-pspy(ib)*sin(vx(ib)*deg2rad)+pspx(ib)*cos(vx(ib)*deg2rad)
1406 
1407  pspxd_tmp=-pspyd(ib)*cos(vx(ib)*deg2rad)-pspxd(ib)*sin(vx(ib)*deg2rad)
1408  pspyd_tmp=-pspyd(ib)*sin(vx(ib)*deg2rad)+pspxd(ib)*cos(vx(ib)*deg2rad)
1409 
1410  fij1=s1(ia,k)+dxa*pspx(ia)+dya*pspy(ia)
1411  fij2=s1(ib,k)+dxb*pspx_tmp+dyb*pspy_tmp
1412 
1413  !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
1414  ! David moved HPRNU and added VHC
1415  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))) !/HPRNU
1416 
1417  txx=0.5_sp*(pspxd(ia)+pspxd_tmp)*viscof
1418  tyy=0.5_sp*(pspyd(ia)+pspyd_tmp)*viscof
1419  ELSE IF(ib == node_northpole)THEN
1420  pspx_tmp=-pspy(ia)*cos(vx(ia)*deg2rad)-pspx(ia)*sin(vx(ia)*deg2rad)
1421  pspy_tmp=-pspy(ia)*sin(vx(ia)*deg2rad)+pspx(ia)*cos(vx(ia)*deg2rad)
1422 
1423  pspxd_tmp=-pspyd(ia)*cos(vx(ia)*deg2rad)-pspxd(ia)*sin(vx(ia)*deg2rad)
1424  pspyd_tmp=-pspyd(ia)*sin(vx(ia)*deg2rad)+pspxd(ia)*cos(vx(ia)*deg2rad)
1425 
1426  fij1=s1(ia,k)+dxa*pspx_tmp+dya*pspy_tmp
1427  fij2=s1(ib,k)+dxb*pspx(ib)+dyb*pspy(ib)
1428 
1429  !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
1430  ! David moved HPRNU and added VHC
1431  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)))
1432 
1433  txx=0.5_sp*(pspxd_tmp+pspxd(ib))*viscof
1434  tyy=0.5_sp*(pspyd_tmp+pspyd(ib))*viscof
1435  END IF
1436 
1437  s1min=minval(s1(nbsn(ia,1:ntsn(ia)-1),k))
1438  s1min=min(s1min, s1(ia,k))
1439  s1max=maxval(s1(nbsn(ia,1:ntsn(ia)-1),k))
1440  s1max=max(s1max, s1(ia,k))
1441  s2min=minval(s1(nbsn(ib,1:ntsn(ib)-1),k))
1442  s2min=min(s2min, s1(ib,k))
1443  s2max=maxval(s1(nbsn(ib,1:ntsn(ib)-1),k))
1444  s2max=max(s2max, s1(ib,k))
1445  IF(fij1 < s1min) fij1=s1min
1446  IF(fij1 > s1max) fij1=s1max
1447  IF(fij2 < s2min) fij2=s2min
1448  IF(fij2 > s2max) fij2=s2max
1449 
1450 ! IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
1451  uij_tmp = -v(i1,k)*cos(xc(i1)*deg2rad)-u(i1,k)*sin(xc(i1)*deg2rad)
1452  vij_tmp = -v(i1,k)*sin(xc(i1)*deg2rad)+u(i1,k)*cos(xc(i1)*deg2rad)
1453 
1454  vx1_tmp = rearth * cos(yije(i,1)*deg2rad) * cos(xije(i,1)*deg2rad)
1455  vy1_tmp = rearth * cos(yije(i,1)*deg2rad) * sin(xije(i,1)*deg2rad)
1456 
1457  vx2_tmp = rearth * cos(yije(i,2)*deg2rad) * cos(xije(i,2)*deg2rad)
1458  vy2_tmp = rearth * cos(yije(i,2)*deg2rad) * sin(xije(i,2)*deg2rad)
1459 
1460  dltxe_tmp = vx2_tmp-vx1_tmp
1461  dltye_tmp = vy2_tmp-vy1_tmp
1462 
1463  fxx=-dtij(i,k)*txx*dltye_tmp
1464  fyy= dtij(i,k)*tyy*dltxe_tmp
1465 
1466  uvn_tmp = vij_tmp*dltxe_tmp - uij_tmp*dltye_tmp
1467  exflux_tmp = -uvn_tmp*dtij(i,k)*((1.0_sp+sign(1.0_sp,uvn_tmp))*fij2+ &
1468  (1.0_sp-sign(1.0_sp,uvn_tmp))*fij1)*0.5_sp
1469 
1470  IF(ia == node_northpole)THEN
1471  xflux(ia,k)=xflux(ia,k)+exflux_tmp+fxx+fyy
1472  xflux_adv(ia,k)=xflux_adv(ia,k)+exflux_tmp
1473  ELSE IF(ib == node_northpole)THEN
1474  xflux(ib,k)=xflux(ib,k)-exflux_tmp-fxx-fyy
1475  xflux_adv(ib,k)=xflux_adv(ib,k)-exflux_tmp
1476  END IF
1477  END IF
1478  END IF
1479  END DO
1480 
1481  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "End: ADV_S_XY(K):",k
1482 
1483  RETURN
1484  END SUBROUTINE adv_s_xy
1485 !==============================================================================|
1486 
1487 !==============================================================================|
1488  SUBROUTINE adv_t_xy(XFLUX,XFLUX_ADV,PTPX,PTPY,PTPXD,PTPYD,VISCOFF,K,CETA)
1489 
1490 !------------------------------------------------------------------------------|
1491 
1492  IMPLICIT NONE
1493  INTEGER, INTENT(IN) :: K
1494  REAL(SP), DIMENSION(0:MT,KB) :: XFLUX,XFLUX_ADV
1495  REAL(SP), DIMENSION(M) :: PTPX,PTPY,PTPXD,PTPYD,VISCOFF
1496  REAL(SP), DIMENSION(3*(NT),KBM1) :: DTIJ
1497  REAL(SP) :: XI,YI
1498  REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2
1499  REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF
1500  REAL(SP) :: FACT,FM1
1501  INTEGER :: I,I1,I2,IA,IB,J,J1,J2,JTMP,JJ,II
1502  REAL(SP) :: TXPI,TYPI
1503 
1504  REAL(SP) :: VX_TMP,VY_TMP,VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP,VX3_TMP,VY3_TMP
1505  REAL(SP) :: XI_TMP,YI_TMP,VXA_TMP,VYA_TMP,VXB_TMP,VYB_TMP
1506  REAL(SP) :: UIJ_TMP,VIJ_TMP,DLTXE_TMP,DLTYE_TMP,UVN_TMP,EXFLUX_TMP
1507  REAL(SP) :: PUPX_TMP,PUPY_TMP,PVPX_TMP,PVPY_TMP
1508  REAL(SP) :: PTPX_TMP,PTPY_TMP,PTPXD_TMP,PTPYD_TMP
1509  REAL(SP) :: U_TMP,V_TMP
1510  REAL(SP) :: X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2
1511 
1512 !! ggao edge calculation
1513  REAL(SP) :: XIJE1_TMP,YIJE1_TMP,XIJE2_TMP,YIJE2_TMP
1514  REAL(SP) :: T1MIN, T1MAX, T2MIN, T2MAX
1515 
1516  REAL(SP) :: CETA
1517 !------------------------------------------------------------------------------!
1518 
1519  ! ALL OTHER PROCESSORS ESCAPE HERE
1520  IF (node_northpole .EQ. 0) RETURN
1521 
1522  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: ADV_T_XY(K):",k
1523 
1524 
1525  SELECT CASE(horizontal_mixing_type)
1526  CASE ('closure')
1527  fact = 1.0_sp
1528  fm1 = 0.0_sp
1529  CASE('constant')
1530  fact = 0.0_sp
1531  fm1 = 1.0_sp
1532  CASE DEFAULT
1533  CALL fatal_error("UNKNOW HORIZONTAL MIXING TYPE:",&
1534  & trim(horizontal_mixing_type) )
1535  END SELECT
1536 
1537 !
1538 !--Initialize Fluxes-----------------------------------------------------------!
1539 !
1540  DO ii=1,npcv
1541  i = ncedge_lst(ii)
1542  ia = niec(i,1)
1543  ib = niec(i,2)
1544  IF(ia == node_northpole)THEN
1545  xflux(ia,k) = 0.0_sp
1546  xflux_adv(ia,k) = 0.0_sp
1547  ELSE IF(ib == node_northpole)THEN
1548  xflux(ib,k) = 0.0_sp
1549  xflux_adv(ib,k) = 0.0_sp
1550  END IF
1551  END DO
1552 
1553 !
1554 !--Loop Over Control Volume Sub-Edges And Calculate Normal Velocity------------!
1555 !
1556  DO ii=1,npcv
1557  i = ncedge_lst(ii)
1558  i1=ntrg(i)
1559  dtij(i,k)=dt1(i1)*dz1(i1,k)
1560  END DO
1561 
1562 !
1563 !--Calculate the Advection and Horizontal Diffusion Terms----------------------!
1564 !
1565  i = node_northpole
1566 
1567  pupx_tmp=0.0_sp
1568  pupy_tmp=0.0_sp
1569  pvpx_tmp=0.0_sp
1570  pvpy_tmp=0.0_sp
1571 
1572  DO j=1,ntve(i)
1573  i1=nbve(i,j)
1574  jtmp=nbvt(i,j)
1575  j1=jtmp+1-(jtmp+1)/4*3
1576  j2=jtmp+2-(jtmp+2)/4*3
1577 
1578  vx_tmp = rearth * cos(vy(i)*deg2rad) * cos(vx(i)*deg2rad) &
1579  * 2._sp /(1._sp+sin(vy(i)*deg2rad))
1580  vy_tmp = rearth * cos(vy(i)*deg2rad) * sin(vx(i)*deg2rad) &
1581  * 2._sp /(1._sp+sin(vy(i)*deg2rad))
1582 
1583  vx1_tmp= rearth * cos(vy(nv(i1,j1))*deg2rad) * cos(vx(nv(i1,j1))*deg2rad) &
1584  * 2._sp /(1._sp+sin(vy(nv(i1,j1))*deg2rad))
1585  vy1_tmp= rearth * cos(vy(nv(i1,j1))*deg2rad) * sin(vx(nv(i1,j1))*deg2rad) &
1586  * 2._sp /(1._sp+sin(vy(nv(i1,j1))*deg2rad))
1587 
1588  vx2_tmp= rearth * cos(yc(i1)*deg2rad) * cos(xc(i1)*deg2rad) &
1589  * 2._sp /(1._sp+sin(yc(i1)*deg2rad))
1590  vy2_tmp= rearth * cos(yc(i1)*deg2rad) * sin(xc(i1)*deg2rad) &
1591  * 2._sp /(1._sp+sin(yc(i1)*deg2rad))
1592 
1593  vx3_tmp= rearth * cos(vy(nv(i1,j2))*deg2rad) * cos(vx(nv(i1,j2))*deg2rad) &
1594  * 2._sp /(1._sp+sin(vy(nv(i1,j2))*deg2rad))
1595  vy3_tmp= rearth * cos(vy(nv(i1,j2))*deg2rad) * sin(vx(nv(i1,j2))*deg2rad) &
1596  * 2._sp /(1._sp+sin(vy(nv(i1,j2))*deg2rad))
1597 
1598  x11=0.5_sp*(vx_tmp+vx1_tmp)
1599  y11=0.5_sp*(vy_tmp+vy1_tmp)
1600  x22=vx2_tmp
1601  y22=vx2_tmp
1602  x33=0.5_sp*(vx_tmp+vx3_tmp)
1603  y33=0.5_sp*(vy_tmp+vy3_tmp)
1604 
1605  u_tmp = -v(i1,k)*cos(xc(i1)*deg2rad)-u(i1,k)*sin(xc(i1)*deg2rad)
1606  v_tmp = -v(i1,k)*sin(xc(i1)*deg2rad)+u(i1,k)*cos(xc(i1)*deg2rad)
1607 
1608  pupx_tmp=pupx_tmp+u_tmp*(y11-y33)
1609  pupy_tmp=pupy_tmp+u_tmp*(x33-x11)
1610  pvpx_tmp=pvpx_tmp+v_tmp*(y11-y33)
1611  pvpy_tmp=pvpy_tmp+v_tmp*(x33-x11)
1612  END DO
1613 
1614  pupx_tmp=pupx_tmp/art1(i)
1615  pupy_tmp=pupy_tmp/art1(i)
1616  pvpx_tmp=pvpx_tmp/art1(i)
1617  pvpy_tmp=pvpy_tmp/art1(i)
1618  tmp1=pupx_tmp**2+pvpy_tmp**2
1619  tmp2=0.5_sp*(pupy_tmp+pvpx_tmp)**2
1620  viscoff(i)=sqrt(tmp1+tmp2)*art1(i)
1621 
1622  IF(k == kbm1) THEN
1623  ah_bottom(i) = (fact*viscoff(i) + fm1)*nn_hvc(i)
1624  END IF
1625 
1626  DO ii=1,npcv
1627  i = ncedge_lst(ii)
1628  i1=ntrg(i)
1629  ia=niec(i,1)
1630  ib=niec(i,2)
1631  IF(ia <= m .AND. ib <= m .AND. i1 <= n)THEN
1632 ! XI=0.5_SP*(XIJE(I,1)+XIJE(I,2))
1633 ! YI=0.5_SP*(YIJE(I,1)+YIJE(I,2))
1634 !! ggao edge calculation
1635  xije1_tmp = rearth * cos(yije(i,1)*deg2rad) * cos(xije(i,1)*deg2rad) &
1636  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
1637  yije1_tmp = rearth * cos(yije(i,1)*deg2rad) * sin(xije(i,1)*deg2rad) &
1638  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
1639 
1640  xije2_tmp = rearth * cos(yije(i,2)*deg2rad) * cos(xije(i,2)*deg2rad) &
1641  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
1642  yije2_tmp = rearth * cos(yije(i,2)*deg2rad) * sin(xije(i,2)*deg2rad) &
1643  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
1644  xi_tmp =0.5_sp*(xije1_tmp+xije2_tmp)
1645  yi_tmp =0.5_sp*(yije1_tmp+yije2_tmp)
1646 
1647  IF(ia == node_northpole .OR. ib == node_northpole)THEN
1648 ! XI_TMP = REARTH * COS(YI*DEG2RAD) * COS(XI*DEG2RAD) &
1649 ! * 2._SP /(1._SP+sin(YI*DEG2RAD))
1650 ! YI_TMP = REARTH * COS(YI*DEG2RAD) * SIN(XI*DEG2RAD) &
1651 ! * 2._SP /(1._SP+sin(YI*DEG2RAD))
1652 
1653  vxa_tmp = rearth * cos(vy(ia)*deg2rad) * cos(vx(ia)*deg2rad) &
1654  * 2._sp /(1._sp+sin(vy(ia)*deg2rad))
1655  vya_tmp = rearth * cos(vy(ia)*deg2rad) * sin(vx(ia)*deg2rad) &
1656  * 2._sp /(1._sp+sin(vy(ia)*deg2rad))
1657 
1658  vxb_tmp = rearth * cos(vy(ib)*deg2rad) * cos(vx(ib)*deg2rad) &
1659  * 2._sp /(1._sp+sin(vy(ib)*deg2rad))
1660  vyb_tmp = rearth * cos(vy(ib)*deg2rad) * sin(vx(ib)*deg2rad) &
1661  * 2._sp /(1._sp+sin(vy(ib)*deg2rad))
1662 
1663 ! IF(IA == NODE_NORTHPOLE)THEN
1664  dxa=xi_tmp-vxa_tmp
1665  dya=yi_tmp-vya_tmp
1666 ! ELSE IF(IB == NODE_NORTHPOLE)THEN
1667  dxb=xi_tmp-vxb_tmp
1668  dyb=yi_tmp-vyb_tmp
1669 ! END IF
1670 ! END IF
1671 
1672  IF(ia == node_northpole)THEN
1673  ptpx_tmp=-ptpy(ib)*cos(vx(ib)*deg2rad)-ptpx(ib)*sin(vx(ib)*deg2rad)
1674  ptpy_tmp=-ptpy(ib)*sin(vx(ib)*deg2rad)+ptpx(ib)*cos(vx(ib)*deg2rad)
1675 
1676  ptpxd_tmp=-ptpyd(ib)*cos(vx(ib)*deg2rad)-ptpxd(ib)*sin(vx(ib)*deg2rad)
1677  ptpyd_tmp=-ptpyd(ib)*sin(vx(ib)*deg2rad)+ptpxd(ib)*cos(vx(ib)*deg2rad)
1678 
1679  fij1=t1(ia,k)+dxa*ptpx(ia)+dya*ptpy(ia)
1680  fij2=t1(ib,k)+dxb*ptpx_tmp+dyb*ptpy_tmp
1681 
1682  !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
1683  ! David moved HPRNU and added VHC
1684  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))) !/HPRNU
1685 
1686  txx=0.5_sp*(ptpxd(ia)+ptpxd_tmp)*viscof
1687  tyy=0.5_sp*(ptpyd(ia)+ptpyd_tmp)*viscof
1688  ELSE IF(ib == node_northpole)THEN
1689  ptpx_tmp=-ptpy(ia)*cos(vx(ia)*deg2rad)-ptpx(ia)*sin(vx(ia)*deg2rad)
1690  ptpy_tmp=-ptpy(ia)*sin(vx(ia)*deg2rad)+ptpx(ia)*cos(vx(ia)*deg2rad)
1691 
1692  ptpxd_tmp=-ptpyd(ia)*cos(vx(ia)*deg2rad)-ptpxd(ia)*sin(vx(ia)*deg2rad)
1693  ptpyd_tmp=-ptpyd(ia)*sin(vx(ia)*deg2rad)+ptpxd(ia)*cos(vx(ia)*deg2rad)
1694 
1695  fij1=t1(ia,k)+dxa*ptpx_tmp+dya*ptpy_tmp
1696  fij2=t1(ib,k)+dxb*ptpx(ib)+dyb*ptpy(ib)
1697 
1698  !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
1699  ! David moved HPRNU and added VHC
1700  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)))
1701 
1702  txx=0.5_sp*(ptpxd_tmp+ptpxd(ib))*viscof
1703  tyy=0.5_sp*(ptpyd_tmp+ptpyd(ib))*viscof
1704  END IF
1705 
1706  t1min=minval(t1(nbsn(ia,1:ntsn(ia)-1),k))
1707  t1min=min(t1min, t1(ia,k))
1708  t1max=maxval(t1(nbsn(ia,1:ntsn(ia)-1),k))
1709  t1max=max(t1max, t1(ia,k))
1710  t2min=minval(t1(nbsn(ib,1:ntsn(ib)-1),k))
1711  t2min=min(t2min, t1(ib,k))
1712  t2max=maxval(t1(nbsn(ib,1:ntsn(ib)-1),k))
1713  t2max=max(t2max, t1(ib,k))
1714  IF(fij1 < t1min) fij1=t1min
1715  IF(fij1 > t1max) fij1=t1max
1716  IF(fij2 < t2min) fij2=t2min
1717  IF(fij2 > t2max) fij2=t2max
1718 
1719 ! IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
1720  uij_tmp = -v(i1,k)*cos(xc(i1)*deg2rad)-u(i1,k)*sin(xc(i1)*deg2rad)
1721  vij_tmp = -v(i1,k)*sin(xc(i1)*deg2rad)+u(i1,k)*cos(xc(i1)*deg2rad)
1722 
1723  vx1_tmp = rearth * cos(yije(i,1)*deg2rad) * cos(xije(i,1)*deg2rad)
1724  vy1_tmp = rearth * cos(yije(i,1)*deg2rad) * sin(xije(i,1)*deg2rad)
1725 
1726  vx2_tmp = rearth * cos(yije(i,2)*deg2rad) * cos(xije(i,2)*deg2rad)
1727  vy2_tmp = rearth * cos(yije(i,2)*deg2rad) * sin(xije(i,2)*deg2rad)
1728 
1729  dltxe_tmp = vx2_tmp-vx1_tmp
1730  dltye_tmp = vy2_tmp-vy1_tmp
1731 
1732  fxx=-dtij(i,k)*txx*dltye_tmp
1733  fyy= dtij(i,k)*tyy*dltxe_tmp
1734 
1735  uvn_tmp = vij_tmp*dltxe_tmp - uij_tmp*dltye_tmp
1736  exflux_tmp = -uvn_tmp*dtij(i,k)*((1.0_sp+sign(1.0_sp,uvn_tmp))*fij2+ &
1737  (1.0_sp-sign(1.0_sp,uvn_tmp))*fij1)*0.5_sp
1738 
1739  IF(ia == node_northpole)THEN
1740  xflux(ia,k)=xflux(ia,k)+exflux_tmp+fxx+fyy
1741  xflux_adv(ia,k)=xflux_adv(ia,k)+exflux_tmp
1742  ELSE IF(ib == node_northpole)THEN
1743  xflux(ib,k)=xflux(ib,k)-exflux_tmp-fxx-fyy
1744  xflux_adv(ib,k)=xflux_adv(ib,k)-exflux_tmp
1745  END IF
1746  END IF
1747  END IF
1748  END DO
1749 
1750  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "End: ADV_T_XY(K):",k
1751 
1752  RETURN
1753  END SUBROUTINE adv_t_xy
1754 !==============================================================================|
1755 
1756  SUBROUTINE adv_n_xy(XFLUX,PWPX,PWPY,ISS,ID,DEP2,SPCDIR,N32)
1757 ! This subroutine is for wave advection at the north pole yzhang
1758 !------------------------------------------------------------------------------|
1759  USE swcomm3, ONLY :mdc,msc
1760 
1761  IMPLICIT NONE
1762 
1763  SAVE
1764 
1765  REAL :: N32(MDC,MSC,0:MT),N32_TMP(MDC,MSC,0:MT)
1766  INTEGER :: ID,ISS,IG,IG2,IDT,IDD ! LWU IG2 IS FOR PLBC
1767  REAL :: CANX,CANY,CANX_TMP,CANY_TMP,ADDEXFLUX2455
1768  REAL :: SPCDIR(MDC,6)
1769  REAL(SP) :: DEP2(MT)
1770  REAL :: DEPLOC,KWAVELOC,CGLOC,NN,ND,SPCSIGL
1771  REAL(SP), DIMENSION(0:MT) :: XFLUX
1772  REAL(SP), DIMENSION(M) :: PWPX,PWPY
1773  REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2
1774  INTEGER :: I,I1,I2,IA,IB,J,J1,J2,JTMP,JJ,II,L
1775  REAL(SP) :: VX_TMP,VY_TMP,VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP,VX3_TMP,VY3_TMP
1776  REAL(SP) :: XI_TMP,YI_TMP,VXA_TMP,VYA_TMP,VXB_TMP,VYB_TMP
1777  REAL(SP) :: UL_DEGREE,DL_DEGREE,CENTER_DEGREE,FF11
1778  REAL(SP) :: UIJ_TMP,VIJ_TMP,DLTXE_TMP,DLTYE_TMP,UVN_TMP,EXFLUX_TMP
1779  REAL(SP) :: PUPX_TMP,PUPY_TMP,PVPX_TMP,PVPY_TMP
1780  REAL(SP) :: PWPX_TMP,PWPY_TMP
1781  REAL(SP) :: U_TMP,V_TMP
1782  REAL(SP) :: ADDYIJE1,ADDYIJE2
1783  REAL(SP) :: ADD_DLTXE ,ADD_DLTYE,ADD_DLTXTRIE,ADD_DLTYTRIE
1784  REAL(SP) :: X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2
1785  REAL(SP) :: XIJE1_TMP,YIJE1_TMP,XIJE2_TMP,YIJE2_TMP
1786  REAL(SP) :: AC1MIN, AC1MAX, AC2MIN, AC2MAX
1787 
1788 !------------------------------------------------------------------------------!
1789  IF (node_northpole .EQ. 0) RETURN
1790 
1791  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: ADV_N_XY"
1792 
1793 ! CAX = 0.0
1794 ! CAY = 0.0
1795 
1796 !
1797 !--Initialize Fluxes-----------------------------------------------------------!
1798 !
1799  DO ii=1,npcv
1800  i = ncedge_lst(ii)
1801  ia = niec(i,1)
1802  ib = niec(i,2)
1803  IF(ia == node_northpole)THEN
1804  xflux(ia) = 0.0_sp
1805  ELSE IF(ib == node_northpole)THEN
1806  xflux(ib) = 0.0_sp
1807  END IF
1808  END DO
1809 !==========YZHANG_NORTHPOLE_TESTING================================
1810  i=node_northpole
1811  pwpx(i)=0.0_sp
1812  pwpy(i)=0.0_sp
1813  add_dltxtrie=0.0_sp
1814  add_dltytrie=0.0_sp
1815  DO j=1,ntsn(i)-1
1816  i1=nbsn(i,j)
1817  i2=nbsn(i,j+1)
1818  ul_degree=22.5
1819  center_degree=0
1820  dl_degree=337.5
1821 
1822  IF (vx(i1)<ul_degree .OR. vx(i1)>dl_degree)THEN
1823  idt=mdc-(mdc*2/8)
1824  idd=id+idt
1825 ! IDD=ID
1826  IF(idd>mdc)THEN
1827  idd=idd-mdc
1828  END IF
1829  n32_tmp(id,iss,i1)=n32(idd,iss,i1)
1830  END IF
1831 
1832  IF (vx(i2)<ul_degree .OR. vx(i2)>dl_degree)THEN
1833  idt=mdc-(mdc*2/8)
1834  idd=id+idt
1835 ! IDD=ID
1836  IF(idd>mdc)THEN
1837  idd=idd-mdc
1838  END IF
1839  n32_tmp(id,iss,i2)=n32(idd,iss,i2)
1840  END IF
1841  DO l=1,7
1842  center_degree=l*45.0
1843  ul_degree=center_degree+22.5
1844  dl_degree=center_degree-22.5
1845 
1846  IF(vx(i1)<ul_degree .AND. vx(i1)>dl_degree) THEN
1847  idt=mdc-(mdc*(l+2)/8)
1848 ! IDT=MDC-(MDC*L/8)
1849  IF(idt<0)THEN
1850  idt=idt+mdc
1851  END IF
1852  idd=id+idt
1853  IF(idd>mdc)THEN
1854  idd=idd-mdc
1855  END IF
1856  n32_tmp(id,iss,i1)=n32(idd,iss,i1)
1857  END IF
1858  IF(vx(i2)<ul_degree .AND. vx(i2)>dl_degree) THEN
1859  idt=mdc-(mdc*(l+2)/8)
1860 ! IDT=MDC-(MDC*L/8)
1861  IF(idt<0)THEN
1862  idt=idt+mdc
1863  END IF
1864  idd=id+idt
1865  IF(idd>mdc)THEN
1866  idd=idd-mdc
1867  END IF
1868  n32_tmp(id,iss,i2)=n32(idd,iss,i2)
1869  END IF
1870  END DO
1871  ff11=0.5*(n32_tmp(id,iss,i1)+n32_tmp(id,iss,i2))
1872  pwpx(i)=pwpx(i)+ff11*dltytrie(i,j)
1873  pwpy(i)=pwpy(i)+ff11*dltxtrie(i,j)
1874  add_dltxtrie=add_dltxtrie+dltxtrie(i,j)
1875  add_dltytrie=add_dltytrie+dltytrie(i,j)
1876  END DO
1877 
1878  pwpx(i)=pwpx(i)/art2(i)
1879  pwpy(i)=pwpy(i)/art2(i)
1880 !============================================================================
1881  addyije1=0.0_sp
1882  addyije2=0.0_sp
1883  DO ii=1,npcv
1884  i = ncedge_lst(ii)
1885  addyije1=addyije1+yije(i,1)
1886  addyije2=addyije2+yije(i,2)
1887  END DO
1888  addyije1=addyije1/npcv
1889  addyije2=addyije2/npcv
1890 
1891  add_dltxe=0.0_sp
1892  add_dltye=0.0_sp
1893 
1894  DO ii=1,npcv
1895  i = ncedge_lst(ii)
1896  i1=ntrg(i)
1897  ia=niec(i,1)
1898  ib=niec(i,2)
1899  IF(ia <= m .AND. ib <= m .AND. i1 <= n)THEN
1900 ! XI=0.5_SP*(XIJE(I,1)+XIJE(I,2))
1901 ! YI=0.5_SP*(YIJE(I,1)+YIJE(I,2))
1902 !! ggao edge calculation
1903  xije1_tmp = rearth * cos(yije(i,1)*deg2rad) * cos(xije(i,1)*deg2rad) &
1904  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
1905  yije1_tmp = rearth * cos(yije(i,1)*deg2rad) * sin(xije(i,1)*deg2rad) &
1906  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
1907 
1908  xije2_tmp = rearth * cos(yije(i,2)*deg2rad) * cos(xije(i,2)*deg2rad) &
1909  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
1910  yije2_tmp = rearth * cos(yije(i,2)*deg2rad) * sin(xije(i,2)*deg2rad) &
1911  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
1912  xi_tmp =0.5_sp*(xije1_tmp+xije2_tmp)
1913  yi_tmp =0.5_sp*(yije1_tmp+yije2_tmp)
1914 
1915  IF(ia == node_northpole .OR. ib == node_northpole)THEN
1916  vxa_tmp = rearth * cos(vy(ia)*deg2rad) * cos(vx(ia)*deg2rad) &
1917  * 2._sp /(1._sp+sin(vy(ia)*deg2rad))
1918  vya_tmp = rearth * cos(vy(ia)*deg2rad) * sin(vx(ia)*deg2rad) &
1919  * 2._sp /(1._sp+sin(vy(ia)*deg2rad))
1920 
1921  vxb_tmp = rearth * cos(vy(ib)*deg2rad) * cos(vx(ib)*deg2rad) &
1922  * 2._sp /(1._sp+sin(vy(ib)*deg2rad))
1923  vyb_tmp = rearth * cos(vy(ib)*deg2rad) * sin(vx(ib)*deg2rad) &
1924  * 2._sp /(1._sp+sin(vy(ib)*deg2rad))
1925 
1926 ! IF(IA == NODE_NORTHPOLE)THEN
1927  dxa=xi_tmp-vxa_tmp
1928  dya=yi_tmp-vya_tmp
1929 ! ELSE IF(IB == NODE_NORTHPOLE)THEN
1930  dxb=xi_tmp-vxb_tmp
1931  dyb=yi_tmp-vyb_tmp
1932 
1933  IF(ia == node_northpole)THEN
1934  pwpx_tmp=-pwpy(ib)*cos(vx(ib)*deg2rad)-pwpx(ib)*sin(vx(ib)*deg2rad)
1935  pwpy_tmp=-pwpy(ib)*sin(vx(ib)*deg2rad)+pwpx(ib)*cos(vx(ib)*deg2rad)
1936 
1937 ! PSPXD_TMP=-PSPYD(IB)*COS(VX(IB)*DEG2RAD)-PSPXD(IB)*SIN(VX(IB)*DEG2RAD)
1938 ! PSPYD_TMP=-PSPYD(IB)*SIN(VX(IB)*DEG2RAD)+PSPXD(IB)*COS(VX(IB)*DEG2RAD)
1939 !======================zhangyang======start=======================
1940  IF (vx(ib)<ul_degree .OR. vx(ib)>dl_degree)THEN
1941  idt=mdc-(mdc*2/8)
1942  idd=id+idt
1943 ! IDD=ID
1944  IF(idd>mdc)THEN
1945  idd=idd-mdc
1946  END IF
1947  n32_tmp(id,iss,ib)=n32(idd,iss,ib)
1948  END IF
1949 
1950  DO l=1,7
1951  center_degree=l*45.0
1952  ul_degree=center_degree+22.5
1953  dl_degree=center_degree-22.5
1954  IF(vx(ib)<ul_degree .AND. vx(ib)>dl_degree) THEN
1955  idt=mdc-(mdc*(l+2)/8)
1956 ! IDT=MDC-(MDC*L/8)
1957  IF(idt<0)THEN
1958  idt=idt+mdc
1959  END IF
1960  idd=id+idt
1961  IF(idd>mdc)THEN
1962  idd=idd-mdc
1963  END IF
1964  n32_tmp(id,iss,ib)=n32(idd,iss,ib)
1965  END IF
1966  END DO
1967  fij1=n32(id,iss,ia)!+DXA*PWPX(IA)+DYA*PWPY(IA)
1968  fij2=n32_tmp(id,iss,ib)!+DXB*PWPX_TMP+DYB*PWPY_TMP
1969 !================zhangyang======end==================================
1970 
1971  !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
1972  ! David moved HPRNU and added VHC
1973 ! VISCOF=(FACT*0.5_SP*(VISCOFF(IA)*NN_HVC(IA)+VISCOFF(IB)*NN_HVC(IB)) +
1974 ! FM1*0.5_SP*(NN_HVC(IA)+NN_HVC(IB))) !/HPRNU
1975 
1976 ! TXX=0.5_SP*(PSPXD(IA)+PSPXD_TMP)*VISCOF
1977 ! TYY=0.5_SP*(PSPYD(IA)+PSPYD_TMP)*VISCOF
1978 
1979  ELSE IF(ib == node_northpole)THEN
1980  pwpx_tmp=-pwpy(ia)*cos(vx(ia)*deg2rad)-pwpx(ia)*sin(vx(ia)*deg2rad)
1981  pwpy_tmp=-pwpy(ia)*sin(vx(ia)*deg2rad)+pwpx(ia)*cos(vx(ia)*deg2rad)
1982 
1983 ! PSPXD_TMP=-PSPYD(IA)*COS(VX(IA)*DEG2RAD)-PSPXD(IA)*SIN(VX(IA)*DEG2RAD)
1984 ! PSPYD_TMP=-PSPYD(IA)*SIN(VX(IA)*DEG2RAD)+PSPXD(IA)*COS(VX(IA)*DEG2RAD)
1985 
1986 !======================zhangyang======start=======================
1987  IF (vx(ia)<ul_degree .OR. vx(ia)>dl_degree)THEN
1988  idt=mdc-(mdc*2/8)
1989  idd=id+idt
1990 ! IDD=ID
1991  IF(idd>mdc)THEN
1992  idd=idd-mdc
1993  END IF
1994  n32_tmp(id,iss,ia)=n32(idd,iss,ia)
1995  END IF
1996 
1997  DO l=1,7
1998  center_degree=l*45.0
1999  ul_degree=center_degree+22.5
2000  dl_degree=center_degree-22.5
2001  IF(vx(ia)<ul_degree .AND. vx(ia)>dl_degree) THEN
2002  idt=mdc-(mdc*(l+2)/8)
2003 ! IDT=MDC-(MDC*L/8)
2004  IF(idt<0)THEN
2005  idt=idt+mdc
2006  END IF
2007  idd=id+idt
2008  IF(idd>mdc)THEN
2009  idd=idd-mdc
2010  END IF
2011  n32_tmp(id,iss,ia)=n32(idd,iss,ia)
2012  END IF
2013  END DO
2014  fij1=n32_tmp(id,iss,ia)!+DXA*PWPX_TMP+DYA*PWPY_TMP
2015  fij2=n32(id,iss,ib)!+DXB*PWPX(IB)+DYB*PWPY(IB)
2016 
2017 !================zhangyang======end==================================
2018  END IF
2019  CALL swapar1(i1,iss,id,dep2(1),kwaveloc,cgloc)
2020  uij_tmp = ua(i1)
2021  vij_tmp = va(i1)
2022 
2023  DO l=1,8
2024  center_degree=l*45.0-22.5
2025  ul_degree=center_degree+22.5
2026  dl_degree=center_degree-22.5
2027 
2028  IF(xc(i1)<ul_degree .AND. xc(i1)>dl_degree) THEN
2029  idt=mdc-(mdc*(l+2)/8-5)
2030 ! IDT=MDC-(MDC*L/8)
2031  IF(idt<0)THEN
2032  idt=idt+mdc
2033  END IF
2034  idd=id+idt
2035  IF(idd>mdc)THEN
2036  idd=idd-mdc
2037  END IF
2038  END IF
2039  END DO
2040 
2041  CALL sproxy(i1,iss,idd,canx,cany,cgloc,spcdir(idd,2),spcdir(idd,3),uij_tmp,vij_tmp)
2042 
2043  vx1_tmp = rearth * cos(yije(i,1)*deg2rad) * cos(xije(i,1)*deg2rad) &
2044  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
2045  vy1_tmp = rearth * cos(yije(i,1)*deg2rad) * sin(xije(i,1)*deg2rad) &
2046  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
2047 
2048  vx2_tmp = rearth * cos(yije(i,2)*deg2rad) * cos(xije(i,2)*deg2rad) &
2049  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
2050  vy2_tmp = rearth * cos(yije(i,2)*deg2rad) * sin(xije(i,2)*deg2rad) &
2051  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
2052 
2053  dltxe_tmp = vx2_tmp-vx1_tmp
2054  dltye_tmp = vy2_tmp-vy1_tmp
2055 
2056  canx_tmp = -cany*cos(xc(i1)*deg2rad)-canx*sin(xc(i1)*deg2rad)
2057  cany_tmp = -cany*sin(xc(i1)*deg2rad)+canx*cos(xc(i1)*deg2rad)
2058 
2059  uvn_tmp = cany_tmp*dltxe_tmp - canx_tmp*dltye_tmp
2060  exflux_tmp = -uvn_tmp*((1.0_sp+sign(1.0_sp,uvn_tmp))*fij2+ &
2061  (1.0_sp-sign(1.0_sp,uvn_tmp))*fij1)*0.5_sp
2062 
2063  IF(ia == node_northpole)THEN
2064  xflux(ia)=xflux(ia)+exflux_tmp
2065 ! XFLUX_ADV(IA,K)=XFLUX_ADV(IA,K)+EXFLUX_TMP
2066  ELSE IF(ib == node_northpole)THEN
2067  xflux(ib)=xflux(ib)-exflux_tmp
2068 ! XFLUX_ADV(IB,K)=XFLUX_ADV(IB,K)-EXFLUX_TMP
2069  END IF
2070  END IF
2071  END IF
2072  END DO
2073 
2074  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "End: ADV_N_XY(ID,ISS):",id,iss
2075 
2076  RETURN
2077  END SUBROUTINE adv_n_xy
2078 
2079 !==============================================================================|
2080 ! CALCULATE THE BAROCLINIC PRESSURE GRADIENT IN SIGMA COORDINATES |
2081 !==============================================================================|
2082 
2083  SUBROUTINE baropg_xy(DRIJK1,DRIJK2)
2084 
2085 !==============================================================================|
2086  IMPLICIT NONE
2087  REAL(SP) :: DRIJK1(0:N,3,KBM1), DRIJK2(0:N,KBM1)
2088  REAL(SP) :: DIJ,DRHO1,DRHO2
2089  INTEGER :: I,II,K,J,J1,J2,IJK
2090  REAL(SP) :: VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP
2091 !==============================================================================|
2092 
2093  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: baropg_XY"
2094 
2095 
2096  DO ii = 1, np
2097  i=np_lst(ii)
2098  DO k=1,kbm1
2099  drhox(i,k)=0.0_sp
2100  drhoy(i,k)=0.0_sp
2101  END DO
2102  END DO
2103 
2104  DO ii = 1, np
2105  i=np_lst(ii)
2106  DO k=1,kbm1
2107  DO j = 1, 3
2108  j1=j+1-int((j+1)/4)*3
2109  j2=j+2-int((j+2)/4)*3
2110  ijk=nbe(i,j)
2111  dij=0.5_sp*(dt(nv(i,j1))+dt(nv(i,j2)))
2112 
2113  vy1_tmp=rearth*cos(vy(nv(i,j1))*deg2rad)*sin(vx(nv(i,j1))*deg2rad)
2114  vy2_tmp=rearth*cos(vy(nv(i,j2))*deg2rad)*sin(vx(nv(i,j2))*deg2rad)
2115 
2116  drho1=(vy1_tmp-vy2_tmp)*drijk1(i,j,k)*dt1(i)
2117  drho2=(vy1_tmp-vy2_tmp)*dij*drijk2(i,k)
2118 
2119  drhox(i,k)=drhox(i,k)+drho1+drho2
2120 
2121  vx1_tmp=rearth*cos(vy(nv(i,j1))*deg2rad)*cos(vx(nv(i,j1))*deg2rad)
2122  vx2_tmp=rearth*cos(vy(nv(i,j2))*deg2rad)*cos(vx(nv(i,j2))*deg2rad)
2123 
2124  drho1=(vx2_tmp-vx1_tmp)*drijk1(i,j,k)*dt1(i)
2125  drho2=(vx2_tmp-vx1_tmp)*dij*drijk2(i,k)
2126 
2127  drhoy(i,k)=drhoy(i,k)+drho1+drho2
2128 
2129  END DO
2130  END DO
2131  END DO
2132 
2133  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "End: baropg_XY"
2134 
2135  RETURN
2136  END SUBROUTINE baropg_xy
2137 !==============================================================================|
2138 
2139 !==============================================================================|
2140  SUBROUTINE shape_coef_xy
2141 
2142 !----------------------------------------------------------------------!
2143 ! This subrountine is used to calculate the coefficient for a linear !
2144 ! function on the x-y plane, i.e.: !
2145 ! r(x,y;phai)=phai_c+cofa1*x+cofa2*y !
2146 ! innc(i)=0 cells on the boundary !
2147 ! innc(i)=1 cells in the interior !
2148 !----------------------------------------------------------------------!
2149 
2150  USE all_vars
2151  IMPLICIT NONE
2152  REAL(DP) X1,X2,X3,Y1,Y2,Y3,DELT,AI1,AI2,AI3,BI1,BI2,BI3,CI1,CI2,CI3
2153  REAL(DP) DELTX,DELTY,TEMP1,ANG1,ANG2,B1,B2,ANGLE
2154  REAL(DP), ALLOCATABLE :: XC_TMP(:),YC_TMP(:),VX_TMP(:),VY_TMP(:)
2155  INTEGER I,II,J,JJ,J1,J2
2156 !
2157 !---------------interior cells-----------------------------------------!
2158 !
2159  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: shape_coef_xy"
2160 
2161  ! ALL OTHER PROCESSORS ESCAPE HERE
2162  IF (node_northpole .EQ. 0) RETURN
2163 
2164 
2165  ALLOCATE(a1u_xy(n,4)); a1u_xy = 0.0_sp
2166  ALLOCATE(a2u_xy(n,4)); a2u_xy = 0.0_sp
2167  ALLOCATE(aw0_xy(n,3)); aw0_xy = 0.0_sp
2168  ALLOCATE(awx_xy(n,3)); awx_xy = 0.0_sp
2169  ALLOCATE(awy_xy(n,3)); awy_xy = 0.0_sp
2170 
2171  ALLOCATE(xc_tmp(0:nt)); xc_tmp = 0.0_sp
2172  ALLOCATE(yc_tmp(0:nt)); yc_tmp = 0.0_sp
2173  ALLOCATE(vx_tmp(0:mt)); vx_tmp = 0.0_sp
2174  ALLOCATE(vy_tmp(0:mt)); vy_tmp = 0.0_sp
2175 
2176  DO i=1,nt
2177  xc_tmp(i) = rearth * cos(yc(i)*deg2rad) * cos(xc(i)*deg2rad) &
2178  * 2._sp /(1._sp+sin(yc(i)*deg2rad))
2179  yc_tmp(i) = rearth * cos(yc(i)*deg2rad) * sin(xc(i)*deg2rad) &
2180  * 2._sp /(1._sp+sin(yc(i)*deg2rad))
2181  END DO
2182 
2183  DO i=1,mt
2184  vx_tmp(i) = rearth * cos(vy(i)*deg2rad) * cos(vx(i)*deg2rad) &
2185  * 2._sp /(1._sp+sin(vy(i)*deg2rad))
2186  vy_tmp(i) = rearth * cos(vy(i)*deg2rad) * sin(vx(i)*deg2rad) &
2187  * 2._sp /(1._sp+sin(vy(i)*deg2rad))
2188  END DO
2189 
2190  DO i=1,n
2191  IF(isbce(i) == 0)THEN
2192  y1 = yc_tmp(nbe(i,1))-yc_tmp(i)
2193  y2 = yc_tmp(nbe(i,2))-yc_tmp(i)
2194  y3 = yc_tmp(nbe(i,3))-yc_tmp(i)
2195  x1=xc_tmp(nbe(i,1))-xc_tmp(i)
2196  x2=xc_tmp(nbe(i,2))-xc_tmp(i)
2197  x3=xc_tmp(nbe(i,3))-xc_tmp(i)
2198 
2199  x1=x1/1000.0_sp
2200  x2=x2/1000.0_sp
2201  x3=x3/1000.0_sp
2202  y1=y1/1000.0_sp
2203  y2=y2/1000.0_sp
2204  y3=y3/1000.0_sp
2205 
2206  delt=(x1*y2-x2*y1)**2+(x1*y3-x3*y1)**2+(x2*y3-x3*y2)**2
2207  delt=delt*1000.0_sp
2208 
2209  a1u_xy(i,1)=(y1+y2+y3)*(x1*y1+x2*y2+x3*y3)- &
2210  (x1+x2+x3)*(y1**2+y2**2+y3**2)
2211  a1u_xy(i,1)=a1u_xy(i,1)/delt
2212  a1u_xy(i,2)=(y1**2+y2**2+y3**2)*x1-(x1*y1+x2*y2+x3*y3)*y1
2213  a1u_xy(i,2)=a1u_xy(i,2)/delt
2214  a1u_xy(i,3)=(y1**2+y2**2+y3**2)*x2-(x1*y1+x2*y2+x3*y3)*y2
2215  a1u_xy(i,3)=a1u_xy(i,3)/delt
2216  a1u_xy(i,4)=(y1**2+y2**2+y3**2)*x3-(x1*y1+x2*y2+x3*y3)*y3
2217  a1u_xy(i,4)=a1u_xy(i,4)/delt
2218 
2219  a2u_xy(i,1)=(x1+x2+x3)*(x1*y1+x2*y2+x3*y3)- &
2220  (y1+y2+y3)*(x1**2+x2**2+x3**2)
2221  a2u_xy(i,1)=a2u_xy(i,1)/delt
2222  a2u_xy(i,2)=(x1**2+x2**2+x3**2)*y1-(x1*y1+x2*y2+x3*y3)*x1
2223  a2u_xy(i,2)=a2u_xy(i,2)/delt
2224  a2u_xy(i,3)=(x1**2+x2**2+x3**2)*y2-(x1*y1+x2*y2+x3*y3)*x2
2225  a2u_xy(i,3)=a2u_xy(i,3)/delt
2226  a2u_xy(i,4)=(x1**2+x2**2+x3**2)*y3-(x1*y1+x2*y2+x3*y3)*x3
2227  a2u_xy(i,4)=a2u_xy(i,4)/delt
2228  end if
2229 
2230  x1=vx_tmp(nv(i,1))-xc_tmp(i)
2231  x2=vx_tmp(nv(i,2))-xc_tmp(i)
2232  x3=vx_tmp(nv(i,3))-xc_tmp(i)
2233  y1=vy_tmp(nv(i,1))-yc_tmp(i)
2234  y2=vy_tmp(nv(i,2))-yc_tmp(i)
2235  y3=vy_tmp(nv(i,3))-yc_tmp(i)
2236 
2237 
2238  ai1=y2-y3
2239  ai2=y3-y1
2240  ai3=y1-y2
2241  bi1=x3-x2
2242  bi2=x1-x3
2243  bi3=x2-x1
2244  ci1=x2*y3-x3*y2
2245  ci2=x3*y1-x1*y3
2246  ci3=x1*y2-x2*y1
2247 
2248  aw0_xy(i,1)=-ci1/2./art(i)
2249  aw0_xy(i,2)=-ci2/2./art(i)
2250  aw0_xy(i,3)=-ci3/2./art(i)
2251  awx_xy(i,1)=-ai1/2./art(i)
2252  awx_xy(i,2)=-ai2/2./art(i)
2253  awx_xy(i,3)=-ai3/2./art(i)
2254  awy_xy(i,1)=-bi1/2./art(i)
2255  awy_xy(i,2)=-bi2/2./art(i)
2256  awy_xy(i,3)=-bi3/2./art(i)
2257  end do
2258 
2259  DEALLOCATE(xc_tmp,yc_tmp,vx_tmp,vy_tmp)
2260 
2261  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "End: shape_coef_xy"
2262 
2263  return
2264  end subroutine shape_coef_xy
2265 
2266 !==============================================================================|
2267 
2268 !==============================================================================|
2269 !==============================================================================|
2270  SUBROUTINE adv_q_xy(XFLUX,PQPX,PQPY,PQPXD,PQPYD,VISCOFF,Q,UQ,VQ,K,UQ1,VQ1,CETA)
2271 
2272 !==============================================================================|
2273 ! Calculate the Turbulent Kinetic Energy and Mixing Length Based on |
2274 ! The Mellor-Yamada Level 2.5 Turbulent Closure Model |
2275 !==============================================================================|
2276 
2277 
2278 !------------------------------------------------------------------------------|
2279 
2280  IMPLICIT NONE
2281  INTEGER, INTENT(IN) :: K
2282  REAL(SP), DIMENSION(0:MT,KB) :: XFLUX,Q
2283  REAL(SP), DIMENSION(M) :: PQPX,PQPY,PQPXD,PQPYD,VISCOFF
2284  REAL(SP), DIMENSION(3*(NT),KBM1) :: DTIJ
2285  REAL(SP) :: XI,YI
2286  REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2
2287  REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF
2288  REAL(SP) :: FACT,FM1
2289  INTEGER :: I,I1,I2,IA,IB,J,J1,J2,JTMP,JJ,II
2290  REAL(SP) :: TXPI,TYPI
2291 
2292  REAL(SP) :: VX_TMP,VY_TMP,VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP,VX3_TMP,VY3_TMP
2293  REAL(SP) :: XI_TMP,YI_TMP,VXA_TMP,VYA_TMP,VXB_TMP,VYB_TMP
2294  REAL(SP) :: UIJ_TMP,VIJ_TMP,DLTXE_TMP,DLTYE_TMP,UVN_TMP,EXFLUX_TMP
2295  REAL(SP) :: PUPX_TMP,PUPY_TMP,PVPX_TMP,PVPY_TMP
2296  REAL(SP) :: PQPX_TMP,PQPY_TMP,PQPXD_TMP,PQPYD_TMP
2297  REAL(SP) :: U_TMP,V_TMP
2298  REAL(SP) :: X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2
2299 
2300  REAL(SP) :: XIJE1_TMP,YIJE1_TMP,XIJE2_TMP,YIJE2_TMP
2301  REAL(SP) :: Q1MIN, Q1MAX, Q2MIN, Q2MAX
2302 
2303  REAL(SP), DIMENSION(0:NT,KB) :: UQ,VQ
2304 
2305  REAL(SP), DIMENSION(0:,:) :: UQ1, VQ1
2306  REAL(SP) :: CETA
2307 !------------------------------------------------------------------------------!
2308  IF (node_northpole .EQ. 0) RETURN
2309 
2310  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "Start: adv_q_xy(K):",k
2311 
2312 
2313  SELECT CASE(horizontal_mixing_type)
2314  CASE ('closure')
2315  fact = 1.0_sp
2316  fm1 = 0.0_sp
2317  CASE('constant')
2318  fact = 0.0_sp
2319  fm1 = 1.0_sp
2320  CASE DEFAULT
2321  CALL fatal_error("UNKNOW HORIZONTAL MIXING TYPE:",&
2322  & trim(horizontal_mixing_type) )
2323  END SELECT
2324 
2325 
2326 !
2327 !--Initialize Fluxes-----------------------------------------------------------!
2328 !
2329  DO ii=1,npcv
2330  i = ncedge_lst(ii)
2331  ia = niec(i,1)
2332  ib = niec(i,2)
2333  IF(ia == node_northpole)THEN
2334  xflux(ia,k) = 0.0_sp
2335  ELSE IF(ib == node_northpole)THEN
2336  xflux(ib,k) = 0.0_sp
2337  END IF
2338  END DO
2339 
2340 !
2341 !--Loop Over Control Volume Sub-Edges And Calculate Normal Velocity------------!
2342 !
2343  DO ii=1,npcv
2344  i = ncedge_lst(ii)
2345  i1=ntrg(i)
2346  dtij(i,k)=dt1(i1)*dz1(i1,k)
2347  END DO
2348 
2349 !
2350 !--Calculate the Advection and Horizontal Diffusion Terms----------------------!
2351 !
2352  i = node_northpole
2353 
2354  IF(i==0) RETURN
2355 
2356  IF(i /= 0)THEN
2357  pupx_tmp=0.0_sp
2358  pupy_tmp=0.0_sp
2359  pvpx_tmp=0.0_sp
2360  pvpy_tmp=0.0_sp
2361 
2362  DO j=1,ntve(i)
2363  i1=nbve(i,j)
2364  jtmp=nbvt(i,j)
2365  j1=jtmp+1-(jtmp+1)/4*3
2366  j2=jtmp+2-(jtmp+2)/4*3
2367 
2368  vx_tmp = rearth * cos(vy(i)*deg2rad) * cos(vx(i)*deg2rad) &
2369  * 2._sp /(1._sp+sin(vy(i)*deg2rad))
2370  vy_tmp = rearth * cos(vy(i)*deg2rad) * sin(vx(i)*deg2rad) &
2371  * 2._sp /(1._sp+sin(vy(i)*deg2rad))
2372 
2373  vx1_tmp= rearth * cos(vy(nv(i1,j1))*deg2rad) * cos(vx(nv(i1,j1))*deg2rad) &
2374  * 2._sp /(1._sp+sin(vy(nv(i1,j1))*deg2rad))
2375  vy1_tmp= rearth * cos(vy(nv(i1,j1))*deg2rad) * sin(vx(nv(i1,j1))*deg2rad) &
2376  * 2._sp /(1._sp+sin(vy(nv(i1,j1))*deg2rad))
2377 
2378  vx2_tmp= rearth * cos(yc(i1)*deg2rad) * cos(xc(i1)*deg2rad) &
2379  * 2._sp /(1._sp+sin(yc(i1)*deg2rad))
2380  vy2_tmp= rearth * cos(yc(i1)*deg2rad) * sin(xc(i1)*deg2rad) &
2381  * 2._sp /(1._sp+sin(yc(i1)*deg2rad))
2382 
2383  vx3_tmp= rearth * cos(vy(nv(i1,j2))*deg2rad) * cos(vx(nv(i1,j2))*deg2rad) &
2384  * 2._sp /(1._sp+sin(vy(nv(i1,j2))*deg2rad))
2385  vy3_tmp= rearth * cos(vy(nv(i1,j2))*deg2rad) * sin(vx(nv(i1,j2))*deg2rad) &
2386  * 2._sp /(1._sp+sin(vy(nv(i1,j2))*deg2rad))
2387 
2388  x11=0.5_sp*(vx_tmp+vx1_tmp)
2389  y11=0.5_sp*(vy_tmp+vy1_tmp)
2390  x22=vx2_tmp
2391  y22=vx2_tmp
2392  x33=0.5_sp*(vx_tmp+vx3_tmp)
2393  y33=0.5_sp*(vy_tmp+vy3_tmp)
2394 
2395  u_tmp = -vq(i1,k)*cos(xc(i1)*deg2rad)-uq(i1,k)*sin(xc(i1)*deg2rad)
2396  v_tmp = -vq(i1,k)*sin(xc(i1)*deg2rad)+uq(i1,k)*cos(xc(i1)*deg2rad)
2397 
2398  pupx_tmp=pupx_tmp+u_tmp*(y11-y33)
2399  pupy_tmp=pupy_tmp+u_tmp*(x33-x11)
2400  pvpx_tmp=pvpx_tmp+v_tmp*(y11-y33)
2401  pvpy_tmp=pvpy_tmp+v_tmp*(x33-x11)
2402  END DO
2403 
2404  pupx_tmp=pupx_tmp/art1(i)
2405  pupy_tmp=pupy_tmp/art1(i)
2406  pvpx_tmp=pvpx_tmp/art1(i)
2407  pvpy_tmp=pvpy_tmp/art1(i)
2408  tmp1=pupx_tmp**2+pvpy_tmp**2
2409  tmp2=0.5_sp*(pupy_tmp+pvpx_tmp)**2
2410  viscoff(i)=sqrt(tmp1+tmp2)*art1(i)
2411 
2412  IF(k == kbm1) THEN
2413  ah_bottom(i) = (fact*viscoff(i) + fm1)*nn_hvc(i)
2414  END IF
2415  endif
2416  DO ii=1,npcv
2417  i = ncedge_lst(ii)
2418  i1=ntrg(i)
2419  ia=niec(i,1)
2420  ib=niec(i,2)
2421 
2422  IF((ia <= m .AND. ib <= m) .AND. i1 <= n)THEN
2423  xije1_tmp = rearth * cos(yije(i,1)*deg2rad) * cos(xije(i,1)*deg2rad) &
2424  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
2425  yije1_tmp = rearth * cos(yije(i,1)*deg2rad) * sin(xije(i,1)*deg2rad) &
2426  * 2._sp /(1._sp+sin(yije(i,1)*deg2rad))
2427 
2428  xije2_tmp = rearth * cos(yije(i,2)*deg2rad) * cos(xije(i,2)*deg2rad) &
2429  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
2430  yije2_tmp = rearth * cos(yije(i,2)*deg2rad) * sin(xije(i,2)*deg2rad) &
2431  * 2._sp /(1._sp+sin(yije(i,2)*deg2rad))
2432  xi_tmp =0.5_sp*(xije1_tmp+xije2_tmp)
2433  yi_tmp =0.5_sp*(yije1_tmp+yije2_tmp)
2434 
2435  IF(ia == node_northpole .OR. ib == node_northpole)THEN
2436  vxa_tmp = rearth * cos(vy(ia)*deg2rad) * cos(vx(ia)*deg2rad) &
2437  * 2._sp /(1._sp+sin(vy(ia)*deg2rad))
2438  vya_tmp = rearth * cos(vy(ia)*deg2rad) * sin(vx(ia)*deg2rad) &
2439  * 2._sp /(1._sp+sin(vy(ia)*deg2rad))
2440 
2441  vxb_tmp = rearth * cos(vy(ib)*deg2rad) * cos(vx(ib)*deg2rad) &
2442  * 2._sp /(1._sp+sin(vy(ib)*deg2rad))
2443  vyb_tmp = rearth * cos(vy(ib)*deg2rad) * sin(vx(ib)*deg2rad) &
2444  * 2._sp /(1._sp+sin(vy(ib)*deg2rad))
2445 
2446  dxa=xi_tmp-vxa_tmp
2447  dya=yi_tmp-vya_tmp
2448  dxb=xi_tmp-vxb_tmp
2449  dyb=yi_tmp-vyb_tmp
2450 
2451  IF(ia == node_northpole)THEN
2452  pqpx_tmp=-pqpy(ib)*cos(vx(ib)*deg2rad)-pqpx(ib)*sin(vx(ib)*deg2rad)
2453  pqpy_tmp=-pqpy(ib)*sin(vx(ib)*deg2rad)+pqpx(ib)*cos(vx(ib)*deg2rad)
2454 
2455  pqpxd_tmp=-pqpyd(ib)*cos(vx(ib)*deg2rad)-pqpxd(ib)*sin(vx(ib)*deg2rad)
2456  pqpyd_tmp=-pqpyd(ib)*sin(vx(ib)*deg2rad)+pqpxd(ib)*cos(vx(ib)*deg2rad)
2457 
2458  fij1=q(ia,k)+dxa*pqpx(ia)+dya*pqpy(ia)
2459  fij2=q(ib,k)+dxb*pqpx_tmp+dyb*pqpy_tmp
2460 
2461  ! VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
2462  ! David moved HPRNU and added VHC
2463  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)))/hprnu
2464 
2465  txx=0.5_sp*(pqpxd(ia)+pqpxd_tmp)*viscof
2466  tyy=0.5_sp*(pqpyd(ia)+pqpyd_tmp)*viscof
2467  ELSE IF(ib == node_northpole)THEN
2468  pqpx_tmp=-pqpy(ia)*cos(vx(ia)*deg2rad)-pqpx(ia)*sin(vx(ia)*deg2rad)
2469  pqpy_tmp=-pqpy(ia)*sin(vx(ia)*deg2rad)+pqpx(ia)*cos(vx(ia)*deg2rad)
2470 
2471  pqpxd_tmp=-pqpyd(ia)*cos(vx(ia)*deg2rad)-pqpxd(ia)*sin(vx(ia)*deg2rad)
2472  pqpyd_tmp=-pqpyd(ia)*sin(vx(ia)*deg2rad)+pqpxd(ia)*cos(vx(ia)*deg2rad)
2473 
2474  fij1=q(ia,k)+dxa*pqpx_tmp+dya*pqpy_tmp
2475  fij2=q(ib,k)+dxb*pqpx(ib)+dyb*pqpy(ib)
2476 
2477  !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
2478  ! David moved HPRNU and added VHC
2479  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)))/hprnu
2480 
2481  txx=0.5_sp*(pqpxd_tmp+pqpxd(ib))*viscof
2482  tyy=0.5_sp*(pqpyd_tmp+pqpyd(ib))*viscof
2483  END IF
2484 
2485  q1min=minval(q(nbsn(ia,1:ntsn(ia)-1),k))
2486  q1min=min(q1min, q(ia,k))
2487  q1max=maxval(q(nbsn(ia,1:ntsn(ia)-1),k))
2488  q1max=max(q1max, q(ia,k))
2489  q2min=minval(q(nbsn(ib,1:ntsn(ib)-1),k))
2490  q2min=min(q2min, q(ib,k))
2491  q2max=maxval(q(nbsn(ib,1:ntsn(ib)-1),k))
2492  q2max=max(q2max, q(ib,k))
2493  IF(fij1 < q1min) fij1=q1min
2494  IF(fij1 > q1max) fij1=q1max
2495  IF(fij2 < q2min) fij2=q2min
2496  IF(fij2 > q2max) fij2=q2max
2497 
2498 ! FXX=-DTIJ(I,K)*TXX*DLTYE(I)
2499 ! FYY= DTIJ(I,K)*TYY*DLTXE(I)
2500 
2501 ! IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
2502  uij_tmp = -vq(i1,k)*cos(xc(i1)*deg2rad)-uq(i1,k)*sin(xc(i1)*deg2rad)
2503  vij_tmp = -vq(i1,k)*sin(xc(i1)*deg2rad)+uq(i1,k)*cos(xc(i1)*deg2rad)
2504 
2505  vx1_tmp = rearth * cos(yije(i,1)*deg2rad) * cos(xije(i,1)*deg2rad)
2506  vy1_tmp = rearth * cos(yije(i,1)*deg2rad) * sin(xije(i,1)*deg2rad)
2507 
2508  vx2_tmp = rearth * cos(yije(i,2)*deg2rad) * cos(xije(i,2)*deg2rad)
2509  vy2_tmp = rearth * cos(yije(i,2)*deg2rad) * sin(xije(i,2)*deg2rad)
2510 
2511  dltxe_tmp = vx2_tmp-vx1_tmp
2512  dltye_tmp = vy2_tmp-vy1_tmp
2513 
2514  fxx=-dtij(i,k)*txx*dltye_tmp
2515  fyy= dtij(i,k)*tyy*dltxe_tmp
2516 
2517  uvn_tmp = vij_tmp*dltxe_tmp - uij_tmp*dltye_tmp
2518  exflux_tmp = -uvn_tmp*dtij(i,k)*((1.0_sp+sign(1.0_sp,uvn_tmp))*fij2+ &
2519  (1.0_sp-sign(1.0_sp,uvn_tmp))*fij1)*0.5_sp
2520 
2521  IF(ia == node_northpole)THEN
2522  xflux(ia,k)=xflux(ia,k)+exflux_tmp+fxx+fyy
2523  ELSE IF(ib == node_northpole)THEN
2524  xflux(ib,k)=xflux(ib,k)-exflux_tmp-fxx-fyy
2525  END IF
2526  END IF
2527  END IF
2528  END DO
2529 
2530  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "End: ADV_Q_XY(K):",k
2531 
2532  RETURN
2533  END SUBROUTINE adv_q_xy
2534 !==============================================================================|
2535 
2536 END MODULE mod_northpole
all_vars::ienode
integer, dimension(:,:), allocatable, target ienode
Definition: mod_main.f90:1029
all_vars::ntsn
integer, dimension(:), allocatable, target ntsn
Definition: mod_main.f90:1023
all_vars::epor
real(sp), dimension(:), allocatable, target epor
Definition: mod_main.f90:1056
all_vars::elrk1
real(sp), dimension(:), allocatable, target elrk1
Definition: mod_main.f90:1121
all_vars::va
real(sp), dimension(:), allocatable, target va
Definition: mod_main.f90:1104
all_vars::d
real(sp), dimension(:), allocatable, target d
Definition: mod_main.f90:1132
all_vars::yije
real(sp), dimension(:,:), allocatable, target yije
Definition: mod_main.f90:1048
mod_northpole::npe
integer npe
Definition: mod_northpole.f90:49
mod_northpole::advave_edge_xy
subroutine advave_edge_xy(XFLUX, YFLUX, IFCETA)
Definition: mod_northpole.f90:202
all_vars::cor
real(sp), dimension(:), allocatable, target cor
Definition: mod_main.f90:1113
all_vars::d1
real(sp), dimension(:), allocatable, target d1
Definition: mod_main.f90:1116
all_vars::adx2d
real(sp), dimension(:), allocatable, target adx2d
Definition: mod_main.f90:1258
mod_northpole::awy_xy
real(dp), dimension(:,:), allocatable awy_xy
Definition: mod_northpole.f90:56
all_vars::psty
real(sp), dimension(:), allocatable, target psty
Definition: mod_main.f90:1255
mod_utils
Definition: mod_utils.f90:40
all_vars::art
real(sp), dimension(:), allocatable, target art
Definition: mod_main.f90:1009
mod_northpole::mp
integer mp
Definition: mod_northpole.f90:49
mod_northpole::find_northpole
subroutine find_northpole
Definition: mod_northpole.f90:62
mod_spherical
Definition: mod_spherical.f90:45
mod_northpole::adv_s_xy
subroutine adv_s_xy(XFLUX, XFLUX_ADV, PSPX, PSPY, PSPXD, PSPYD, VISCOFF, K, CETA)
Definition: mod_northpole.f90:1215
all_vars::uark
real(sp), dimension(:), allocatable, target uark
Definition: mod_main.f90:1108
all_vars::el
real(sp), dimension(:), allocatable, target el
Definition: mod_main.f90:1134
all_vars::advua
real(sp), dimension(:), allocatable, target advua
Definition: mod_main.f90:1256
all_vars::v
real(sp), dimension(:,:), allocatable, target v
Definition: mod_main.f90:1269
mod_northpole::baropg_xy
subroutine baropg_xy(DRIJK1, DRIJK2)
Definition: mod_northpole.f90:2084
mod_utils::dbg_set
logical function dbg_set(vrb)
Definition: mod_utils.f90:182
mod_northpole::npedge_lst
integer, dimension(:), allocatable npedge_lst
Definition: mod_northpole.f90:52
mod_northpole::node_northpole
integer node_northpole
Definition: mod_northpole.f90:48
mod_northpole::awx_xy
real(dp), dimension(:,:), allocatable awx_xy
Definition: mod_northpole.f90:56
all_vars::art1
real(sp), dimension(:), allocatable, target art1
Definition: mod_main.f90:1010
all_vars::yc
real(sp), dimension(:), allocatable, target yc
Definition: mod_main.f90:1004
all_vars::t1
real(sp), dimension(:,:), allocatable, target t1
Definition: mod_main.f90:1307
mod_par::ngid_x
integer, dimension(:), pointer ngid_x
Definition: mod_par.f90:63
all_vars::cc_hvc
real(sp), dimension(:), allocatable cc_hvc
Definition: mod_main.f90:1302
all_vars::w
real(sp), dimension(:,:), allocatable, target w
Definition: mod_main.f90:1279
mod_obcs::iucp
real(sp), dimension(:), allocatable iucp
Definition: mod_obcs.f90:110
all_vars::egf
real(sp), dimension(:), allocatable, target egf
Definition: mod_main.f90:1136
all_vars::dltytrie
real(sp), dimension(:,:), allocatable, target dltytrie
Definition: mod_main.f90:1063
all_vars::xije
real(sp), dimension(:,:), allocatable, target xije
Definition: mod_main.f90:1047
mod_par
Definition: mod_par.f90:40
mod_northpole::adv_t_xy
subroutine adv_t_xy(XFLUX, XFLUX_ADV, PTPX, PTPY, PTPXD, PTPYD, VISCOFF, K, CETA)
Definition: mod_northpole.f90:1489
all_vars::art2
real(sp), dimension(:), allocatable, target art2
Definition: mod_main.f90:1011
mod_northpole::vertvl_edge_xy
subroutine vertvl_edge_xy(XFLUX, CETA)
Definition: mod_northpole.f90:1134
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::isbc
integer, dimension(:), allocatable, target isbc
Definition: mod_main.f90:1026
swcomm3::mdc
integer mdc
Definition: swmod1.f90:1672
all_vars::s1
real(sp), dimension(:,:), allocatable, target s1
Definition: mod_main.f90:1308
all_vars::iec
integer, dimension(:,:), allocatable, target iec
Definition: mod_main.f90:1028
mod_northpole
Definition: mod_northpole.f90:40
mod_northpole::adv_n_xy
subroutine adv_n_xy(XFLUX, PWPX, PWPY, ISS, ID, DEP2, SPCDIR, N32)
Definition: mod_northpole.f90:1757
all_vars::drhox
real(sp), dimension(:,:), allocatable, target drhox
Definition: mod_main.f90:1326
all_vars::wubot
real(sp), dimension(:), allocatable, target wubot
Definition: mod_main.f90:1185
all_vars::niec
integer, dimension(:,:), allocatable, target niec
Definition: mod_main.f90:1032
all_vars::nbvt
integer, dimension(:,:), allocatable, target nbvt
Definition: mod_main.f90:1036
swapar1
subroutine swapar1(I, IS, ID, DEP2, KWAVEL, CGOL)
Definition: swancom5.f90:92
all_vars::pstx
real(sp), dimension(:), allocatable, target pstx
Definition: mod_main.f90:1254
mod_northpole::cell_northarea
integer, dimension(:), allocatable cell_northarea
Definition: mod_northpole.f90:51
all_vars::vx
real(sp), dimension(:), allocatable, target vx
Definition: mod_main.f90:1001
mod_northpole::aw0_xy
real(dp), dimension(:,:), allocatable aw0_xy
Definition: mod_northpole.f90:56
all_vars::grav_e
real(sp), dimension(:), allocatable, target grav_e
Definition: mod_main.f90:1013
all_vars::wvbot
real(sp), dimension(:), allocatable, target wvbot
Definition: mod_main.f90:1186
all_vars::wvsurf2
real(sp), dimension(:), allocatable, target wvsurf2
Definition: mod_main.f90:1183
all_vars::vaf
real(sp), dimension(:), allocatable, target vaf
Definition: mod_main.f90:1106
all_vars::nn_hvc
real(sp), dimension(:), allocatable nn_hvc
Definition: mod_main.f90:1303
mod_northpole::shape_coef_xy
subroutine shape_coef_xy
Definition: mod_northpole.f90:2141
all_vars::vy
real(sp), dimension(:), allocatable, target vy
Definition: mod_main.f90:1002
all_vars::ntve
integer, dimension(:), allocatable, target ntve
Definition: mod_main.f90:1022
all_vars::nbe
integer, dimension(:,:), allocatable, target nbe
Definition: mod_main.f90:1020
mod_northpole::adv_uv_edge_xy
subroutine adv_uv_edge_xy(XFLUX, YFLUX, CETA, STG, K_STG)
Definition: mod_northpole.f90:671
all_vars::nv
integer, dimension(:,:), allocatable, target nv
Definition: mod_main.f90:1018
all_vars::xijc
real(sp), dimension(:), allocatable, target xijc
Definition: mod_main.f90:1044
mod_northpole::a2u_xy
real(dp), dimension(:,:), allocatable a2u_xy
Definition: mod_northpole.f90:55
mod_northpole::find_cellside
subroutine find_cellside
Definition: mod_northpole.f90:154
mod_northpole::np_lst
integer, dimension(:), allocatable np_lst
Definition: mod_northpole.f90:54
all_vars::drhoy
real(sp), dimension(:,:), allocatable, target drhoy
Definition: mod_main.f90:1327
mod_northpole::adv_q_xy
subroutine adv_q_xy(XFLUX, PQPX, PQPY, PQPXD, PQPYD, VISCOFF, Q, UQ, VQ, K, UQ1, VQ1, CETA)
Definition: mod_northpole.f90:2271
mod_northpole::npcv
integer npcv
Definition: mod_northpole.f90:49
all_vars::wusurf2
real(sp), dimension(:), allocatable, target wusurf2
Definition: mod_main.f90:1182
all_vars::ua
real(sp), dimension(:), allocatable, target ua
Definition: mod_main.f90:1103
mod_northpole::ncedge_lst
integer, dimension(:), allocatable ncedge_lst
Definition: mod_northpole.f90:53
all_vars::dz
real(sp), dimension(:,:), allocatable, target dz
Definition: mod_main.f90:1092
all_vars
Definition: mod_main.f90:959
all_vars::ady2d
real(sp), dimension(:), allocatable, target ady2d
Definition: mod_main.f90:1259
sproxy
subroutine sproxy(I1, IS, ID, CAXL, CAYL, CG0L, ECOSL, ESINL, UX2L, UY2L)
Definition: swancom5.f90:146
mod_northpole::extel_edge_xy
subroutine extel_edge_xy(K, XFLUX)
Definition: mod_northpole.f90:996
all_vars::yijc
real(sp), dimension(:), allocatable, target yijc
Definition: mod_main.f90:1045
mod_northpole::a1u_xy
real(dp), dimension(:,:), allocatable a1u_xy
Definition: mod_northpole.f90:55
all_vars::dltxtrie
real(sp), dimension(:,:), allocatable, target dltxtrie
Definition: mod_main.f90:1064
all_vars::nbve
integer, dimension(:,:), allocatable, target nbve
Definition: mod_main.f90:1034
all_vars::dt1
real(sp), dimension(:), allocatable, target dt1
Definition: mod_main.f90:1117
mod_obcs
Definition: mod_obcs.f90:62
mod_northpole::node_northarea
integer, dimension(:), allocatable node_northarea
Definition: mod_northpole.f90:50
swcomm3::msc
integer msc
Definition: swmod1.f90:1673
mod_utils::fatal_error
subroutine fatal_error(ER1, ER2, ER3, ER4)
Definition: mod_utils.f90:230
all_vars::h1
real(sp), dimension(:), allocatable, target h1
Definition: mod_main.f90:1115
mod_northpole::extuv_edge_xy
subroutine extuv_edge_xy(K)
Definition: mod_northpole.f90:1073
all_vars::xc
real(sp), dimension(:), allocatable, target xc
Definition: mod_main.f90:1003
all_vars::uaf
real(sp), dimension(:), allocatable, target uaf
Definition: mod_main.f90:1105
all_vars::dry2d
real(sp), dimension(:), allocatable, target dry2d
Definition: mod_main.f90:1261
all_vars::dz1
real(sp), dimension(:,:), allocatable, target dz1
Definition: mod_main.f90:1096
all_vars::elf1
real(sp), dimension(:), allocatable, target elf1
Definition: mod_main.f90:1123
all_vars::isbce
integer, dimension(:), allocatable, target isbce
Definition: mod_main.f90:1027
mod_northpole::np
integer np
Definition: mod_northpole.f90:49
all_vars::nbsn
integer, dimension(:,:), allocatable, target nbsn
Definition: mod_main.f90:1030
all_vars::ah_bottom
real(sp), dimension(:), allocatable, target ah_bottom
Definition: mod_main.f90:1345
mod_utils::dbg_sbr
integer, parameter dbg_sbr
Definition: mod_utils.f90:69
all_vars::vark
real(sp), dimension(:), allocatable, target vark
Definition: mod_main.f90:1109
mod_northpole::advection_edge_xy
subroutine advection_edge_xy(XFLUX, YFLUX)
Definition: mod_northpole.f90:459
mod_northpole::mp_lst
integer, dimension(:), allocatable mp_lst
Definition: mod_northpole.f90:54
all_vars::advva
real(sp), dimension(:), allocatable, target advva
Definition: mod_main.f90:1257
all_vars::dt
real(sp), dimension(:), allocatable, target dt
Definition: mod_main.f90:1133
mod_utils::dbg_log
integer, parameter dbg_log
Definition: mod_utils.f90:65
all_vars::drx2d
real(sp), dimension(:), allocatable, target drx2d
Definition: mod_main.f90:1260
Generated by   doxygen 1.8.15