My Project
mod_action_ex.f90
Go to the documentation of this file.
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 MODULE mod_action_ex
13 
14  IMPLICIT NONE
15  SAVE
16 
17  REAL, ALLOCATABLE :: n31(:,:)
18  REAL, ALLOCATABLE :: n32(:,:,:)
19 
20  CONTAINS
21 !
22 !==============================================================================|
23 !
24 !==========================================================================|
25 !
26  SUBROUTINE action_allo
27 
28  USE all_vars, ONLY : mt
29  USE swcomm3, ONLY : mdc,msc
30 ! USE MOD_USGRID, ONLY : MDC,MSC
31  IMPLICIT NONE
32 
33  ALLOCATE(n31(mdc,msc)) ; n31 = 0.0
34  ALLOCATE(n32(mdc,msc,0:mt)); n32 = 0.0
35 
36  RETURN
37  END SUBROUTINE action_allo
38 !
39 !==========================================================================|
40 !
41 !==========================================================================|
42 !
43  SUBROUTINE action_deallo
44 
45  IMPLICIT NONE
46 
47  DEALLOCATE(n31)
48  DEALLOCATE(n32)
49 
50  RETURN
51  END SUBROUTINE action_deallo
52 !
53 !==========================================================================|
54 !
55 !==========================================================================|
56 !
57  SUBROUTINE algorithm_fct(CAS,IG,DTW,IDCMIN,IDCMAX)
58 
59  USE swcomm3, ONLY : mdc,msc
60  USE m_genarr, ONLY : spcsig
61 ! USE MOD_USGRID, ONLY : SPCSIG,MDC,MSC
62  USE vars_wave, ONLY : mt,ac2
63 ! USE ALL_VARS, ONLY : MT,AC2
64  IMPLICIT NONE
65 
66  INTEGER :: ISS,ID,IG
67  REAL :: CASR,CASL,FLUX1,FLUX2,FLUXLP,FLUXLM,FLUXHP,FLUXHM
68  REAL :: MIN11,MIN22,MIN33,ADLP,ADLM
69  REAL, DIMENSION(MDC,MSC) :: ADP,ADM,NL
70  REAL :: CAS(MDC,MSC,10)
71  REAL :: DTW,DS
72 
73  INTEGER, DIMENSION(MSC) :: IDCMIN,IDCMAX
74 
75  n31 = 0.0
76 
77  DO iss = 1,msc
78  IF(iss == 1)THEN
79  ds = spcsig(iss+1) - spcsig(iss)
80  DO id = 1,mdc
81  casr = 0.5*(cas(id,iss,1)+cas(id,iss+1,1))
82  flux1 = 0.5*(casr+abs(casr))*ac2(id,iss,ig)
83  flux2 = 0.5*(casr-abs(casr))*ac2(id,iss+1,ig)
84  fluxlp = flux1+flux2
85 
86  fluxlm = 0.0
87 
88  fluxhp = casr*0.5*(ac2(id,iss,ig)+ac2(id,iss+1,ig))
89  fluxhm = 0.0
90 
91  nl(id,iss) = ac2(id,iss,ig)-(fluxlp-fluxlm)*dtw/ds
92 
93  adp(id,iss) = fluxhp-fluxlp
94  adm(id,iss) = fluxhm-fluxlm
95  END DO
96  ELSE IF(iss == msc)THEN
97  ds = spcsig(iss) - spcsig(iss-1)
98  DO id = 1,mdc
99  casr = cas(id,iss,1)
100  flux1 = 0.5*(casr+abs(casr))*ac2(id,iss,ig)
101  flux2 = 0.0
102  fluxlp = flux1+flux2
103 
104  casl = cas(id,iss-1,1)
105  flux1 = 0.5*(casl+abs(casl))*ac2(id,iss-1,ig)
106  flux2 = 0.5*(casl-abs(casl))*ac2(id,iss,ig)
107  fluxlm = flux1+flux2
108 
109  fluxhp = casr*ac2(id,iss,ig)
110  fluxhm = casl*ac2(id,iss-1,ig)
111 
112  nl(id,iss) = ac2(id,iss,ig)-(fluxlp-fluxlm)*dtw/ds
113 
114  adp(id,iss) = fluxhp-fluxlp
115  adm(id,iss) = fluxhm-fluxlm
116  END DO
117  ELSE
118  ds = spcsig(iss) - spcsig(iss-1)
119  DO id = 1,mdc
120  casr = 0.5*(cas(id,iss,1)+cas(id,iss+1,1))
121  flux1 = 0.5*(casr+abs(casr))*ac2(id,iss,ig)
122  flux2 = 0.5*(casr-abs(casr))*ac2(id,iss+1,ig)
123  fluxlp = flux1+flux2
124 
125  casl = 0.5*(cas(id,iss,1)+cas(id,iss-1,1))
126  flux1 = 0.5*(casl+abs(casl))*ac2(id,iss-1,ig)
127  flux2 = 0.5*(casl-abs(casl))*ac2(id,iss,ig)
128  fluxlm = flux1+flux2
129 
130  fluxhp = casr*0.5*(ac2(id,iss,ig)+ac2(id,iss+1,ig))
131  fluxhm = casl*0.5*(ac2(id,iss,ig)+ac2(id,iss-1,ig))
132 
133  nl(id,iss) = ac2(id,iss,ig)-(fluxlp-fluxlm)*dtw/ds
134 
135  adp(id,iss) = fluxhp-fluxlp
136  adm(id,iss) = fluxhm-fluxlm
137  END DO
138  END IF
139  END DO
140 
141  DO iss = 1,msc
142  IF(iss == 1)THEN
143  ds = spcsig(iss+1) - spcsig(iss)
144  DO id = 1,mdc
145  min11 = abs(adp(id,iss))
146  min22 = sign(1.,adp(id,iss))*(nl(id,iss+2)-nl(id,iss+1))*ds/dtw
147  adlp = min(min11,min22)
148  adlp = max(0.,adlp)
149  adlp = sign(1.,adp(id,iss))*adlp
150 
151  min11 = abs(adm(id,iss))
152  min22 = sign(1.,adm(id,iss))*(nl(id,iss+1)-nl(id,iss))*ds/dtw
153  adlm = min(min11,min22)
154  adlm = max(0.,adlm)
155  adlm = sign(1.,adm(id,iss))*adlm
156 
157  n31(id,iss) = nl(id,iss)-(adlp-adlm)*dtw/ds
158  END DO
159  ELSE IF(iss == 2)THEN
160  ds = spcsig(iss) - spcsig(iss-1)
161  DO id = 1,mdc
162  min11 = abs(adp(id,iss))
163  min22 = sign(1.,adp(id,iss))*(nl(id,iss+2)-nl(id,iss+1))*ds/dtw
164  min33 = sign(1.,adp(id,iss))*(nl(id,iss)-nl(id,iss-1))*ds/dtw
165  adlp = min(min11,min22,min33)
166  adlp = max(0.,adlp)
167  adlp = sign(1.,adp(id,iss))*adlp
168 
169  min11 = abs(adm(id,iss))
170  min22 = sign(1.,adm(id,iss))*(nl(id,iss+1)-nl(id,iss))*ds/dtw
171  adlm = min(min11,min22)
172  adlm = max(0.,adlm)
173  adlm = sign(1.,adm(id,iss))*adlm
174 
175  n31(id,iss) = nl(id,iss)-(adlp-adlm)*dtw/ds
176  END DO
177  ELSE IF(iss == msc-1)THEN
178  ds = spcsig(iss) - spcsig(iss-1)
179  DO id = 1,mdc
180  min11 = abs(adp(id,iss))
181  min33 = sign(1.,adp(id,iss))*(nl(id,iss)-nl(id,iss-1))*ds/dtw
182  adlp = min(min11,min33)
183  adlp = max(0.,adlp)
184  adlp = sign(1.,adp(id,iss))*adlp
185 
186  min11 = abs(adm(id,iss))
187  min22 = sign(1.,adm(id,iss))*(nl(id,iss+1)-nl(id,iss))*ds/dtw
188  min33 = sign(1.,adm(id,iss))*(nl(id,iss-1)-nl(id,iss-2))*ds/dtw
189  adlm = min(min11,min22,min33)
190  adlm = max(0.,adlm)
191  adlm = sign(1.,adm(id,iss))*adlm
192 
193  n31(id,iss) = nl(id,iss)-(adlp-adlm)*dtw/ds
194  END DO
195  ELSE IF(iss == msc)THEN
196  ds = spcsig(iss) - spcsig(iss-1)
197  DO id = 1,mdc
198  min11 = abs(adp(id,iss))
199  min33 = sign(1.,adp(id,iss))*(nl(id,iss)-nl(id,iss-1))*ds/dtw
200  adlp = min(min11,min33)
201  adlp = max(0.,adlp)
202  adlp = sign(1.,adp(id,iss))*adlp
203 
204  min11 = abs(adm(id,iss))
205  min33 = sign(1.,adm(id,iss))*(nl(id,iss-1)-nl(id,iss-2))*ds/dtw
206  adlm = min(min11,min33)
207  adlm = max(0.,adlm)
208  adlm = sign(1.,adm(id,iss))*adlm
209 
210  n31(id,iss) = nl(id,iss)-(adlp-adlm)*dtw/ds
211  END DO
212  ELSE
213  ds = spcsig(iss) - spcsig(iss-1)
214  DO id = 1,mdc
215  min11 = abs(adp(id,iss))
216  min22 = sign(1.,adp(id,iss))*(nl(id,iss+2)-nl(id,iss+1))*ds/dtw
217  min33 = sign(1.,adp(id,iss))*(nl(id,iss)-nl(id,iss-1))*ds/dtw
218  adlp = min(min11,min22,min33)
219  adlp = max(0.,adlp)
220  adlp = sign(1.,adp(id,iss))*adlp
221 
222  min11 = abs(adm(id,iss))
223  min22 = sign(1.,adm(id,iss))*(nl(id,iss+1)-nl(id,iss))*ds/dtw
224  min33 = sign(1.,adm(id,iss))*(nl(id,iss-1)-nl(id,iss-2))*ds/dtw
225  adlm = min(min11,min22,min33)
226  adlm = max(0.,adlm)
227  adlm = sign(1.,adm(id,iss))*adlm
228 
229  n31(id,iss) = nl(id,iss)-(adlp-adlm)*dtw/ds
230  END DO
231  END IF
232  END DO
233 
234  RETURN
235  END SUBROUTINE algorithm_fct
236 !============================================================================|
237 !============================================================================|
238  SUBROUTINE algorithm_crank_nicolson(CAD,IG,DTW,IDCMIN,IDCMAX,DD)
239 
240  USE swcomm3, ONLY : mdc,msc
241 ! USE MOD_USGRID, ONLY : MDC,MSC
242  IMPLICIT NONE
243  INTEGER :: ISS,ID,IDM1,IDM2,IDP1,MDCM,IG,II
244  INTEGER :: IDDUM
245  INTEGER, DIMENSION(MSC) :: IDCMIN,IDCMAX
246  REAL, PARAMETER :: ZETA = 0.5
247  REAL :: CAD(:,:,:)
248  REAL :: DTW,DD
249  REAL :: N32M,N32P
250 
251  REAL,DIMENSION(MDC) :: A,B,C,R,U
252 
253  DO iss = 1,msc
254  DO iddum = idcmin(iss),idcmax(iss)
255  id = mod(iddum-1+mdc,mdc)+1
256  idp1 = mod(iddum+mdc,mdc)+1
257  idm1 = mod(iddum-2+mdc,mdc)+1
258 
259  b(id) = 1.0
260  IF(id == 1)THEN
261  a(id) = 0.0
262  ELSE
263  a(id) = -0.5*zeta*dtw*cad(idm1,iss,1)/dd
264  END IF
265  IF(id == mdc)THEN
266  c(id) = 0.0
267  ELSE
268  c(id) = 0.5*zeta*dtw*cad(idp1,iss,1)/dd
269  END IF
270 
271  IF(id == 1)THEN
272  r(id) = cad(idp1,iss,1)*n31(idp1,iss)
273  r(id) = (1.0-zeta)*0.5*dtw*r(id)/dd
274  r(id) = n31(id,iss)-r(id)
275  ELSE IF(id == mdc)THEN
276  r(id) = -cad(idm1,iss,1)*n31(idm1,iss)
277  r(id) = (1.0-zeta)*0.5*dtw*r(id)/dd
278  r(id) = n31(id,iss)-r(id)
279  ELSE
280  r(id) = cad(idp1,iss,1)*n31(idp1,iss)-cad(idm1,iss,1)*n31(idm1,iss)
281  r(id) = (1.0-zeta)*0.5*dtw*r(id)/dd
282  r(id) = n31(id,iss)-r(id)
283  END IF
284 
285  END DO
286 
287  CALL tridag(a,b,c,r,u,mdc)
288 
289  DO iddum = idcmin(iss),idcmax(iss)
290  id = mod(iddum-1+mdc,mdc)+1
291  n32(id,iss,ig) = u(id)
292  END DO
293  END DO
294 
295  RETURN
296  END SUBROUTINE algorithm_crank_nicolson
297 !==========================================================================|
298 !
299 !==================================================================================!
300  SUBROUTINE tridag(A,B,C,R,U,N)
301  IMPLICIT NONE
302  INTEGER :: N,J
303  REAL,DIMENSION(N) :: A,B,C,R,U
304  INTEGER, PARAMETER :: NMAX = 500
305  REAL BET,GAM(NMAX)
306 
307  IF(b(1) == 0.)pause 'TRIDAG: REWRITE EQUATIONS'
308  bet = b(1)
309  u(1) = r(1)/bet
310  DO j=2,n
311  gam(j) = c(j-1)/bet
312  bet = b(j)-a(j)*gam(j)
313  IF(bet == 0.)pause 'TRIDAG FAILED'
314  u(j) = (r(j)-a(j)*u(j-1))/bet
315  END DO
316  DO j=n-1,1,-1
317  u(j) = u(j)-gam(j+1)*u(j+1)
318  END DO
319 
320  RETURN
321  END SUBROUTINE tridag
322 !==========================================================================|
323 !
324 !==========================================================================|
325 
326  SUBROUTINE adv_n(DTW)
327 
328 !------------------------------------------------------------------------------|
329 
330  USE vars_wave
331 ! USE ALL_VARS
332  USE swcomm3
333  USE m_genarr
334 
335  IMPLICIT NONE
336 
337  REAL(SP),ALLOCATABLE :: EXFLUX(:,:),UN(:,:),PSPX(:,:,:),PSPY(:,:,:),FF1(:,:),FIJ1(:,:,:)
338  REAL(SP),ALLOCATABLE :: XFLUX(:,:,:)
339  REAL(SP),ALLOCATABLE :: XFLUX_TMP(:)
340  REAL(SP) :: UTMP,VTMP,UN1,X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2,XI,YI
341  REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1_TMP
342  REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF,TEMP,STPOINT
343  REAL(SP) :: FACT,FM1
344  INTEGER :: I,I1,I2,IA,IB,J,J1,J2,K,JTMP,JJ,II
345  REAL(SP) :: AC1MIN, AC1MAX, AC2MIN, AC2MAX
346  INTEGER :: ID,ISS,IG,IG2
347  REAL(SP) :: XIN,YIN,XIC,YIC,CANX,CANY
348  REAL(SP) :: CAX(MDC,MSC),CAY(MDC,MSC)
349 ! REAL(SP) :: DTW,RF(MDC,MSC),DF(MDC,MSC)
350  REAL(SP) :: RF(MDC,MSC),DF(MDC,MSC)
351  REAL :: DTW
352  REAL(SP) :: SPCDIR2(MDC),SPCDIR3(MDC),CANX_2D(MDC,MSC),CANY_2D(MDC,MSC),CGO_1D(MSC),CGO_2D(MDC,MSC)
353  INTEGER :: NTSN_T
354  REAL(SP) :: N32_T,XFLUX_T
355  REAL(SP),ALLOCATABLE :: DEP2(:),AC2LOC(:)
356  REAL(SP) :: DEPLOC,KWAVELOC,CGLOC,NN,ND,SPCSIGL
357  REAL(SP) :: UA_NODE,VA_NODE
358  INTEGER :: CNT
359 
360  REAL,ALLOCATABLE ::EXFLUXA(:,:),EXFLUXB(:,:),UNA(:,:),UNBB(:,:),FIJ2(:,:,:),&
361  FIJ1A(:,:,:),FIJ1B(:,:,:),FIJ2A(:,:,:),FIJ2B(:,:,:), &
362  N32_TMPP2(:,:),N32_TMPP3(:,:)
363  REAL :: DLTXEA,DLTXEB,DLTXETMP,DLTYETMP,C_ICE
364  REAL :: CANXAB(MDC,MSC),CANYA(MDC,MSC),CANYB(MDC,MSC)
365  REAL(SP) :: XFLUX_1D(0:MT),PSPX_1D(M),PSPY_1D(M)
366  REAL :: RATIO(MSC)
367  REAL :: EXFLUX_TMPA(MDC,MSC),EXFLUX_TMPB(MDC,MSC)
368  REAL :: EXFLUX_TMP(MDC,MSC)
369 
370 !------------------------------------------------------------------------------!
371  ALLOCATE(dep2(mt),ac2loc(0:mt),xflux_tmp(0:mt))
372  ALLOCATE(exflux(mdc,msc),ff1(mdc,msc),fij1(mdc,msc,2),un(mdc,msc))
373  ALLOCATE(xflux(mdc,msc,0:mt));xflux=0.0
374  ALLOCATE(pspx(mdc,msc,m));pspx=0.0
375  ALLOCATE(pspy(mdc,msc,m));pspy=0.0
376  IF(high_latitude_wave)THEN
377  ALLOCATE(exfluxa(mdc,msc),exfluxb(mdc,msc),fij2(mdc,msc,0:m),&
378  fij1a(mdc,msc,0:m),fij1b(mdc,msc,0:m),fij2a(mdc,msc,0:m),fij2b(mdc,msc,0:m),una(mdc,msc),unbb(mdc,msc))
379  ALLOCATE(n32_tmpp2(mdc,msc),n32_tmpp3(mdc,msc))
380  END IF
381  dep2(1:mt) = compda(1:mt,jdp2)
382 
383  cax = 0.0
384  cay = 0.0
385 
386 ! DO ISS = 1,MSC
387 ! DO ID = 1,MDC
388 !!
389 !!--Initialize Fluxes-----------------------------------------------------------!
390 !!
391 ! XFLUX = 0.0
392 ! PSPX = 0.0
393 ! PSPY = 0.0
394 !# if defined(PLBC)
395 ! call replace_N32(N32,ID,ISS)
396 !# endif
397 !
398  IF(.NOT. high_latitude_wave)THEN
399 
400  DO i = 1,m
401  DO j=1,ntsn(i)-1
402  i1=nbsn(i,j)
403  i2=nbsn(i,j+1)
404 !!$ IF(DEP2(I1) <= DEPMIN .AND. DEP2(I2) > DEPMIN)THEN
405 !!$ FF1=0.5*(N32(ID,ISS,I)+N32(ID,ISS,I2))
406 !!$ ELSE IF(DEP2(I1) > DEPMIN .AND. DEP2(I2) <= DEPMIN)THEN
407 !!$ FF1=0.5*(N32(ID,ISS,I1)+N32(ID,ISS,I))
408 !!$ ELSE IF(DEP2(I1) <= DEPMIN .AND. DEP2(I2) <= DEPMIN)THEN
409 !!$ FF1=N32(ID,ISS,I)
410 !!$ ELSE
411 !!$ FF1=0.5*(N32(ID,ISS,I1)+N32(ID,ISS,I2))
412 !!$ END IF
413  ff1(:,:)=0.5*(n32(:,:,i1)+n32(:,:,i2))
414 !!$# if defined (SPHERICAL)
415 !!$ XTMP = VX(I2)*TPI-VX(I1)*TPI
416 !!$ XTMP1 = VX(I2)-VX(I1)
417 !!$ IF(XTMP1 > 180.0)THEN
418 !!$ XTMP = -360.0*TPI+XTMP
419 !!$ ELSE IF(XTMP1 < -180.0)THEN
420 !!$ XTMP = 360.0*TPI+XTMP
421 !!$ END IF
422 !!$ TXPI=XTMP*COS(DEG2RAD*VY(I))
423 !!$ TYPI=(VY(I1)-VY(I2))*TPI
424 !!$ PSPX(I)=PSPX(I)+FF1*TYPI
425 !!$ PSPY(I)=PSPY(I)+FF1*TXPI
426 !!$# else
427 !!$ PSPX(I)=PSPX(I)+FF1*(VY(I1)-VY(I2))
428 !!$ PSPY(I)=PSPY(I)+FF1*(VX(I2)-VX(I1))
429 !!$# endif
430  pspx(:,:,i)=pspx(:,:,i)+ff1(:,:)*dltytrie(i,j)
431  pspy(:,:,i)=pspy(:,:,i)+ff1(:,:)*dltxtrie(i,j)
432  END DO
433  pspx(:,:,i)=pspx(:,:,i)/art2(i)
434  pspy(:,:,i)=pspy(:,:,i)/art2(i)
435  END DO
436 !# if defined (PLBC)
437 ! PSPY=0.0_SP
438 !# endif
439  DO i=1,ncv_i
440  i1 = ntrg(i)
441  ia = niec(i,1)
442  ib = niec(i,2)
443 !!$ XI = 0.5*(XIJE(I,1)+XIJE(I,2))
444 !!$ YI = 0.5*(YIJE(I,1)+YIJE(I,2))
445 !!$
446 !!$# if defined (SPHERICAL)
447 !!$ X1_DP=XIJE(I,1)
448 !!$ Y1_DP=YIJE(I,1)
449 !!$ X2_DP=XIJE(I,2)
450 !!$ Y2_DP=YIJE(I,2)
451 !!$ XII = XCG2(I)
452 !!$ YII = YCG2(I)
453 !!$ XI=XII
454 !!$ XTMP = XI*TPI-VX(IA)*TPI
455 !!$ XTMP1 = XI-VX(IA)
456 !!$ IF(XTMP1 > 180.0)THEN
457 !!$ XTMP = -360.0*TPI+XTMP
458 !!$ ELSE IF(XTMP1 < -180.0)THEN
459 !!$ XTMP = 360.0*TPI+XTMP
460 !!$ END IF
461 !!$
462 !!$ DXA=XTMP*VAL_COS_VY(IA)
463 !!$ DYA=(YI-VY(IA))*TPI
464 !!$ XTMP = XI*TPI-VX(IB)*TPI
465 !!$ XTMP1 = XI-VX(IB)
466 !!$ IF(XTMP1 > 180.0)THEN
467 !!$ XTMP = -360.0*TPI+XTMP
468 !!$ ELSE IF(XTMP1 < -180.0)THEN
469 !!$ XTMP = 360.0*TPI+XTMP
470 !!$ END IF
471 !!$
472 !!$ DXB=XTMP*VAL_COS_VY(IB)
473 !!$ DYB=(YI-VY(IB))*TPI
474 !!$# else
475 !!$ DXA = XI - VX(IA)
476 !!$ DYA = YI - VY(IA)
477 !!$ DXB = XI - VX(IB)
478 !!$ DYB = YI - VY(IB)
479 !!$# endif
480 !!$
481 !!$ FIJ1=N32(ID,ISS,IA)+DXA*PSPX(IA)+DYA*PSPY(IA)
482 !!$ FIJ2=N32(ID,ISS,IB)+DXB*PSPX(IB)+DYB*PSPY(IB)
483 
484 ! DEVELOPMENT TESTING - FIRST ORDER WAVE ACTION ADVECTION IN GEOG. SPACE
485  fij1(:,:,1)=n32(:,:,ia)+dltxncve(i,1)*pspx(:,:,ia)+dltyncve(i,1)*pspy(:,:,ia)
486  fij1(:,:,2)=n32(:,:,ib)+dltxncve(i,2)*pspx(:,:,ib)+dltyncve(i,2)*pspy(:,:,ib)
487 
488 !DO ISS=1,MSC
489 !DO ID=1,MDC
490 ! AC1MIN=MINVAL(N32(ID,ISS,NBSN(IA,1:NTSN(IA)-1)))
491 ! AC1MIN=MIN(AC1MIN, N32(ID,ISS,IA))
492 ! AC1MAX=MAXVAL(N32(ID,ISS,NBSN(IA,1:NTSN(IA)-1)))
493 ! AC1MAX=MAX(AC1MAX, N32(ID,ISS,IA))
494 ! AC2MIN=MINVAL(N32(ID,ISS,NBSN(IB,1:NTSN(IB)-1)))
495 ! AC2MIN=MIN(AC2MIN, N32(ID,ISS,IB))
496 ! AC2MAX=MAXVAL(N32(ID,ISS,NBSN(IB,1:NTSN(IB)-1)))
497 ! AC2MAX=MAX(AC2MAX, N32(ID,ISS,IB))
498 ! IF(FIJ1(ID,ISS,1) < AC1MIN) FIJ1(ID,ISS,1) = AC1MIN
499 ! IF(FIJ1(ID,ISS,1) > AC1MAX) FIJ1(ID,ISS,1) = AC1MAX
500 ! IF(FIJ1(ID,ISS,2) < AC2MIN) FIJ1(ID,ISS,2) = AC2MIN
501 ! IF(FIJ1(ID,ISS,2) > AC2MAX) FIJ1(ID,ISS,2) = AC2MAX
502 !END DO
503 !END DO
504 
505  deploc = (dep2(nv(i1,1))+dep2(nv(i1,2))+dep2(nv(i1,3)))/3.0
506  IF(deploc <= depmin)THEN
507 ! *** depth is negative ***
508 ! KWAVEL = -1.
509  cgo_1d = 0.
510  ELSE
511 ! *** call KSCIP1 to compute KWAVE and CGO ***
512  CALL kscip2(spcsig,deploc,cgo_1d)
513  END IF
514 
515 ! CALL SWAPAR1(I1,ISS,ID,DEP2(1),KWAVELOC,CGLOC)
516 
517  utmp = (compda(nv(i1,1),jvx2)+compda(nv(i1,2),jvx2)+ &
518  compda(nv(i1,3),jvx2))/3.0
519  vtmp = (compda(nv(i1,1),jvy2)+compda(nv(i1,2),jvy2)+ &
520  compda(nv(i1,3),jvy2))/3.0
521 
522  DO id=1,mdc
523  cgo_2d(id,:)=cgo_1d
524  END DO
525  spcdir2=spcdir(:,2)
526  spcdir3=spcdir(:,3)
527  CALL sproxy2(canx_2d ,cany_2d , &
528  cgo_2d ,spcdir2,spcdir3,utmp ,vtmp )
529 
530  un = cany_2d*dltxe(i) - canx_2d*dltye(i)
531 !# if defined (PLBC)
532 ! UN = - CANX*DLTYE(I)
533 !# endif
534  exflux=-un*((1.0+sign(1.0,un))*fij1(:,:,2)+(1.0-sign(1.0,un))*fij1(:,:,1))*0.5
535 
536  xflux(:,:,ia)=xflux(:,:,ia)+exflux
537  xflux(:,:,ib)=xflux(:,:,ib)-exflux
538 
539  END DO
540 
541  ELSE
542 
543 
544  END IF
545 
546  DO iss=1,msc
547  DO id=1,mdc
548  DO i = 1,m
549  IF(isonb_w(i) /= 0)THEN
550 !# if !defined(PLBC)
551  deploc = dep2(i)
552 !# else
553 ! CALL replace_node_wave(I,IG2)
554 ! !print*,'I,IG2===',I,IG2
555 ! DEPLOC = DEP2(IG2)
556 !# endif
557  IF(deploc <= depmin)THEN
558 ! *** depth is negative ***
559  kwaveloc = -1.
560  cgloc = 0.
561  ELSE
562 ! *** call KSCIP1 to compute KWAVE and CGO ***
563  spcsigl = spcsig(iss)
564  CALL kscip1(1,spcsigl,deploc,kwaveloc,cgloc,nn,nd)
565  ENDIF
566  cax(id,iss) = cgloc * spcdir(id,2)
567  cay(id,iss) = cgloc * spcdir(id,3)
568 !
569 ! --- adapt the velocities in case of diffraction
570 !
571  IF(idiffr == 1 .AND. pdiffr(3) /= 0.)THEN
572 !JQI CAX(ID,ISS) = CAX(ID,ISS)*DIFPARAM(I)
573 !JQI CAY(ID,ISS) = CAY(ID,ISS)*DIFPARAM(I)
574  END IF
575 !
576 ! --- ambient currents added
577 !
578  IF(icur == 1)THEN
579  cax(id,iss) = cax(id,iss) + compda(i,jvx2)
580  cay(id,iss) = cay(id,iss) + compda(i,jvy2)
581  END IF
582 
583  canx = cax(id,iss)
584  cany = cay(id,iss)
585 ! write(100,*)I,ID,ISS,CANX,CANY
586 
587  fij1_tmp = n32(id,iss,i)
588 
589  IF(nbsn(i,ntsn(i)-1) > m)THEN
590 
591  un1 = cany*(vx(i)-vx(nbsn(i,2))) &
592  -canx*(vy(i)-vy(nbsn(i,2)))
593 !# if defined (PLBC)
594 ! UN1 =-CANX*(VY(I)-VY(NBSN(I,2)))
595 !# endif
596  ELSE IF(nbsn(i,2) > m)THEN
597  un1 = cany*(vx(nbsn(i,ntsn(i)-1))-vx(i)) &
598  -canx*(vy(nbsn(i,ntsn(i)-1))-vy(i))
599 !# if defined (PLBC)
600 ! UN1 = -CANX*(VY(NBSN(I,NTSN(I)-1))-VY(I))
601 !# endif
602  ELSE
603  un1 = cany*(vx(nbsn(i,ntsn(i)-1))-vx(nbsn(i,2))) &
604  -canx*(vy(nbsn(i,ntsn(i)-1))-vy(nbsn(i,2)))
605 !# if defined (PLBC)
606 ! UN1 =-CANX*(VY(NBSN(I,NTSN(I)-1))-VY(NBSN(I,2)))
607 !# endif
608  END IF
609 
610  un1 = 0.5*un1
611 
612 ! EXFLUX = MAX(0.0,-UN1*FIJ1_TMP)
613  xflux(id,iss,i) = xflux(id,iss,i)+max(0.0,-un1*fij1_tmp)
614  END IF
615  END DO
616  END DO
617  END DO
618 !
619 !-Accumulate Fluxes at Boundary Nodes
620 !
621  DO iss=1,msc
622  DO id=1,mdc
623  END DO
624  END DO
625 
626  DO i = 1,m
627 !--Update Action Density ------------------------------------------------------!
628 !
629  IF(dep2(i) > depmin)THEN
630  ac2(:,:,i) = n32(:,:,i)-xflux(:,:,i)/art1(i)*dtw
631  ac2(:,:,i) = max(0.0,ac2(:,:,i))
632  ELSE
633  ac2(:,:,i) = 0.0_sp
634  END IF
635  END DO
636 
637 !# if defined(PLBC)
638 ! CALL replace_ac2(ID,ISS)
639 !# endif
640 
641  DO iss=1,msc
642  DO id=1,mdc
643  END DO
644  END DO
645 
646 !----yzhang----
647 !----yzhang----
648  DEALLOCATE(dep2,ac2loc,xflux_tmp,exflux,ff1,fij1,un)
649  DEALLOCATE(xflux,pspx,pspy)
650  IF(high_latitude_wave)THEN
651  DEALLOCATE(exfluxa,exfluxb,fij2)
652  DEALLOCATE(fij1a,fij1b,fij2a,fij2b,una,unbb)
653  DEALLOCATE(n32_tmpp2,n32_tmpp3)
654  END IF
655 
656  RETURN
657  END SUBROUTINE adv_n
658 
659 !==============================================================================|
660  SUBROUTINE adv_hl(FF1,I,I1,I2)
661 ! This subroutine is for wave advection at high latitude,
662 ! in order to solve the invalid scalar assumption problem.
663 ! yzhang
664 !==========YZHANG_NORTHPOLE_TESTING================================
665  USE vars_wave
666  USE swcomm3
667  USE m_genarr
668  USE mod_northpole
669 
670  IMPLICIT NONE
671  REAL(SP) :: FF1(MDC,MSC),N32_TMPP(MDC,MSC),N32_TMPP4(MDC,MSC)
672  INTEGER :: ISS, ID,I1,I2,IDT,IDD,III,I,IDD1,IDD2
673  REAL :: UL_DEGREE,CENTER_DEGREE,DL_DEGREE,DIFLAT1,DIFLAT2,RATELAT1,RATELAT2
674 
675 
676  DO iss = 1,msc
677  DO id = 1,mdc
678  IF(i1==node_northpole)THEN
679  ul_degree=22.5
680  center_degree=0
681  dl_degree=337.5
682 
683  IF (vx(i)<ul_degree .OR. vx(i)>dl_degree)THEN
684  idt=(mdc*2/8)
685  idd=id+idt
686 ! IDD=ID
687  IF(idd>mdc)THEN
688  idd=idd-mdc
689  END IF
690  n32_tmpp(id,iss)=n32(idd,iss,i1)
691  END IF
692  DO iii=1,7
693  center_degree=iii*45.0
694  ul_degree=center_degree+22.5
695  dl_degree=center_degree-22.5
696  IF(vx(i)<ul_degree .AND. vx(i)>dl_degree) THEN
697  idt=(mdc*(iii+2)/8)
698 ! IDT=(MDC*III/8)
699  IF(idt>=mdc)THEN
700  idt=idt-mdc
701  END IF
702  idd=id+idt
703  IF(idd>mdc)THEN
704  idd=idd-mdc
705  END IF
706  n32_tmpp(id,iss)=n32(idd,iss,i1)
707  END IF
708  END DO
709  ff1(id,iss)=0.5*(n32_tmpp(id,iss)+n32(id,iss,i2))
710  ELSE IF(i2==node_northpole)THEN
711  ul_degree=22.5
712  center_degree=0
713  dl_degree=337.5
714  IF (vx(i)<ul_degree .OR. vx(i)>dl_degree)THEN
715  idt=(mdc*2/8)
716  idd=id+idt
717 ! IDD=ID
718  IF(idd>mdc)THEN
719  idd=idd-mdc
720  END IF
721  n32_tmpp(id,iss)=n32(idd,iss,i2)
722  END IF
723  DO iii=1,7
724  center_degree=iii*45.0
725  ul_degree=center_degree+22.5
726  dl_degree=center_degree-22.5
727  IF(vx(i)<ul_degree .AND. vx(i)>dl_degree) THEN
728  idt=(mdc*(iii+2)/8)
729 ! IDT=(MDC*III/8)
730  IF(idt>=mdc)THEN
731  idt=idt-mdc
732  END IF
733  idd=id+idt
734  IF(idd>mdc)THEN
735  idd=idd-mdc
736  END IF
737  n32_tmpp(id,iss)=n32(idd,iss,i2)
738  END IF
739  END DO
740  ff1(id,iss)=0.5*(n32_tmpp(id,iss)+n32(id,iss,i1))
741  ELSE IF(vy(i)>80.AND.i1/=node_northpole.AND.i2/=node_northpole)THEN
742  diflat1=vx(i)-vx(i1)
743  diflat2=vx(i)-vx(i2)
744  IF(diflat1 > 180.0)THEN
745  diflat1 = -360.0+diflat1
746  ELSE IF(diflat1 < -180.0)THEN
747  diflat1 = 360.0+diflat1
748  END IF
749  IF(diflat2 > 180.0)THEN
750  diflat2 = -360.0+diflat2
751  ELSE IF(diflat2 < -180.0)THEN
752  diflat2 = 360.0+diflat2
753  END IF
754  ratelat1=diflat1/360.
755  ratelat2=diflat2/360.
756  idd1=anint(ratelat1*mdc)+id
757  idd2=anint(ratelat2*mdc)+id
758  IF(idd1<=0)THEN
759  idd1=idd1+mdc
760  END IF
761 
762  IF(idd1>mdc)THEN
763  idd1=idd1-mdc
764  END IF
765 
766  IF(idd2<=0)THEN
767  idd2=idd2+mdc
768  END IF
769 
770  IF(idd2>mdc)THEN
771  idd2=idd2-mdc
772  END IF
773  n32_tmpp(id,iss)=n32(idd1,iss,i1)
774  n32_tmpp4(id,iss)=n32(idd2,iss,i2)
775  ff1(id,iss)=0.5*(n32_tmpp(id,iss)+n32_tmpp4(id,iss))
776 ! ELSE
777 ! FF1(ID,ISS,I)=0.5*(N32(ID,ISS,I1)+N32(ID,ISS,I2))
778  END IF
779  END DO
780  END DO
781 
782  RETURN
783  END SUBROUTINE adv_hl
784 
785 !=========================================================================
786  SUBROUTINE adv_hl_fij(N32_TMPP3,N32_TMPP2,I,IA,IB)
787 ! This subroutine is for wave energy trasport at high latitudes
788 ! yzhang
789 !=========================================================================
790  USE vars_wave
791  USE swcomm3
792  USE m_genarr
793  USE mod_northpole
794 
795  IMPLICIT NONE
796  REAL ::FF1(MDC,MSC),N32_TMPP2(MDC,MSC),N32_TMPP3(MDC,MSC)
797  INTEGER :: ISS, ID,IA,IB,IDT,IDD,III,I,IDD1,IDD2
798  REAL :: UL_DEGREE,CENTER_DEGREE,DL_DEGREE,DIFLAT1,DIFLAT2,RATELAT1,RATELAT2
799 
800 
801  n32_tmpp2=n32(:,:,ia)
802  n32_tmpp3=n32(:,:,ib)
803 
804  DO iss=1,msc
805  DO id=1,mdc
806  IF(ia==node_northpole)THEN
807  ul_degree=22.5
808  center_degree=0
809  dl_degree=337.5
810  IF (vx(ib)<ul_degree .OR. vx(ib)>dl_degree)THEN
811  idt=(mdc*2/8)
812  idd=id+idt
813 ! IDD=ID
814  IF(idd>mdc)THEN
815  idd=idd-mdc
816  END IF
817  n32_tmpp2(id,iss)=n32(idd,iss,ia)
818  END IF
819  DO iii=1,7
820  center_degree=iii*45.0
821  ul_degree=center_degree+22.5
822  dl_degree=center_degree-22.5
823  IF(vx(ib)<ul_degree .AND. vx(ib)>dl_degree) THEN
824  idt=(mdc*(iii+2)/8)
825 ! IDT=(MDC*III/8)
826  IF(idt>=mdc)THEN
827  idt=idt-mdc
828  END IF
829  idd=id+idt
830  IF(idd>mdc)THEN
831  idd=idd-mdc
832  END IF
833  n32_tmpp2(id,iss)=n32(idd,iss,ia)
834  END IF
835  END DO
836  ELSE IF(ib==node_northpole)THEN
837  ul_degree=22.5
838  center_degree=0
839  dl_degree=337.5
840  IF (vx(ia)<ul_degree .OR. vx(ia)>dl_degree)THEN
841  idt=(mdc*2/8)
842  idd=id+idt
843 ! IDD=ID
844  IF(idd>mdc)THEN
845  idd=idd-mdc
846  END IF
847  n32_tmpp2(id,iss)=n32(idd,iss,ib)
848  END IF
849  DO iii=1,7
850  center_degree=iii*45.0
851  ul_degree=center_degree+22.5
852  dl_degree=center_degree-22.5
853  IF(vx(ia)<ul_degree .AND. vx(ia)>dl_degree) THEN
854  idt=(mdc*(iii+2)/8)
855 ! IDT=(MDC*III/8)
856  IF(idt>=mdc)THEN
857  idt=idt-mdc
858  END IF
859  idd=id-idt
860  IF(idd<1)THEN
861  idd=idd+mdc
862  END IF
863  n32_tmpp2(id,iss)=n32(idd,iss,ia)
864  END IF
865  END DO
866  ELSE IF((vy(ia)>80.OR.vy(ib)>80).AND.ib/=node_northpole.AND.ia/=node_northpole)THEN
867  diflat1=vx(ib)-vx(ia)
868  IF(diflat1 > 180.0)THEN
869  diflat1 = -360.0+diflat1
870  ELSE IF(diflat1 < -180.0)THEN
871  diflat1 = 360.0+diflat1
872  END IF
873  ratelat1=diflat1/360.
874  idd1=anint(ratelat1*mdc)+id
875  IF(idd1<=0)THEN
876  idd1=idd1+mdc
877  END IF
878 
879  IF(idd1>mdc)THEN
880  idd1=idd1-mdc
881  END IF
882 
883  n32_tmpp2(id,iss)=n32(idd1,iss,ia)
884 ! ELSE
885 ! N32_TMPP2(ID,ISS,IA)=N32(ID,ISS,IA)
886  END IF
887  IF(ib==node_northpole)THEN
888  ul_degree=22.5
889  center_degree=0
890  dl_degree=337.5
891  IF (vx(ia)<ul_degree .OR. vx(ia)>dl_degree)THEN
892  idt=(mdc*2/8)
893  idd=id+idt
894 ! IDD=ID
895  IF(idd>mdc)THEN
896  idd=idd-mdc
897  END IF
898  n32_tmpp3(id,iss)=n32(idd,iss,ib)
899  END IF
900  DO iii=1,7
901  center_degree=iii*45.0
902  ul_degree=center_degree+22.5
903  dl_degree=center_degree-22.5
904  IF(vx(ia)<ul_degree .AND. vx(ia)>dl_degree) THEN
905  idt=(mdc*(iii+2)/8)
906 ! IDT=(MDC*III/8)
907  IF(idt>=mdc)THEN
908  idt=idt-mdc
909  END IF
910  idd=id+idt
911  IF(idd>mdc)THEN
912  idd=idd-mdc
913  END IF
914  n32_tmpp3(id,iss)=n32(idd,iss,ib)
915  END IF
916  END DO
917  ELSE IF(ia==node_northpole)THEN
918  ul_degree=22.5
919  center_degree=0
920  dl_degree=337.5
921  IF (vx(ib)<ul_degree .OR. vx(ib)>dl_degree)THEN
922  idt=(mdc*2/8)
923  idd=id+idt
924 ! IDD=ID
925  IF(idd>mdc)THEN
926  idd=idd-mdc
927  END IF
928  n32_tmpp3(id,iss)=n32(idd,iss,ia)
929  END IF
930  DO iii=1,7
931  center_degree=iii*45.0
932  ul_degree=center_degree+22.5
933  dl_degree=center_degree-22.5
934  IF(vx(ib)<ul_degree .AND. vx(ib)>dl_degree) THEN
935  idt=(mdc*(iii+2)/8)
936 ! IDT=(MDC*III/8)
937  IF(idt>=mdc)THEN
938  idt=idt-mdc
939  END IF
940  idd=id-idt
941  IF(idd<1)THEN
942  idd=idd+mdc
943  END IF
944  n32_tmpp3(id,iss)=n32(idd,iss,ib)
945  END IF
946  END DO
947  ELSE IF((vy(ia)>80.OR.vy(ib)>80).AND.ib/=node_northpole.AND.ia/=node_northpole)THEN
948  diflat1=vx(ia)-vx(ib)
949  IF(diflat1 > 180.0)THEN
950  diflat1 = -360.0+diflat1
951  ELSE IF(diflat1 < -180.0)THEN
952  diflat1 = 360.0+diflat1
953  END IF
954  ratelat1=diflat1/360.
955  idd1=anint(ratelat1*mdc)+id
956  IF(idd1<=0)THEN
957  idd1=idd1+mdc
958  END IF
959 
960  IF(idd1>mdc)THEN
961  idd1=idd1-mdc
962  END IF
963 
964  n32_tmpp3(id,iss)=n32(idd1,iss,ib)
965 ! ELSE
966 ! N32_TMPP3(ID,ISS,IB)=N32(ID,ISS,IB)
967  END IF
968 ! FIJ1A(ID,ISS,IA)=N32(ID,ISS,IA)+DLTXNCVE(I,1)*PSPX(ID,ISS,IA)+DLTYNCVE(I,1)*PSPY(ID,ISS,IA)
969 ! FIJ2A(ID,ISS,IB)=N32_TMPP2(ID,ISS,IB)+DLTXNCVE(I,2)*PSPX(ID,ISS,IB)+DLTYNCVE(I,2)*PSPY(ID,ISS,IB)
970  END DO
971  END DO
972 
973  RETURN
974  END SUBROUTINE adv_hl_fij
975 
976 !==========================================================================
977 !
978  SUBROUTINE actt2(IA,IB,RATIO)
979 ! This subroutine is for Ice induced wave attenuation
980 ! The parameter 'RATIO' is the dimentional wave attenuation coefficient
981 ! yzhang
982 !
983 !****************************************************************
984 !
985  USE m_genarr
986  USE swcomm2
987  USE swcomm3
988  USE swcomm4
989  USE ocpcomm4
990  USE m_parall
991  USE all_vars
992  USE vars_wave
993  USE ice_state
994 
995 ! USE ice_state , only :aice , vice!, floe_diam
996 ! USE ICE_THERM_VERTICAL ,only:hicen_old
997 ! USE ice_model_size ,only: ncat
998 
999  IMPLICIT NONE
1000  REAL :: C_ICE,COES(7,9),COES2(7,9),ICETHI,SPCSIGL,DEPLOC,TWAVE
1001  REAL :: ICEATT,ICEATTS(7),ICEATT1,ICEATT2,RATE,RATIO(MSC)
1002  REAL :: mean_diam,mean_diam_u,mean_diam_d,floe_diam
1003  REAL ,parameter :: ee=2,ff=0.9,dmin=20
1004  REAL(SP):: DEP3(MT),DIS,DTMP_DP,X1_DP,X2_DP,Y1_DP,Y2_DP
1005  INTEGER :: I,ISS,ID,ni,mm,IA,IB
1006 
1007  dep3(1:mt) = compda(1:mt,jdp2)
1008  coes2(1,:)=(/0.0,0.0,0.0,0.0,-0.003509,0.1473,-2.121,11.24,-22.42/)
1009  coes2(2,:)=(/0.0,0.0,0.0,0.0,-0.003723,0.173,-2.833,18.33,-43.89/)
1010  coes2(3,:)=(/0.0,0.0,0.0,0.0,-0.001264,0.07156,-1.364,9.616,-25.34/)
1011  coes2(4,:)=(/0.0,0.0,0.0,0.0,0.001401,-0.05012,0.5915,-3.329,5.372/)
1012  coes2(5,:)=(/0.0,0.0,0.0,0.0,0.0009987,-0.0403,0.5441,-3.471,6.658/)
1013  coes2(6,:)=(/0.0,0.0,0.0,0.0,0.0002636,-0.01267,0.1853,-1.403,2.459/)
1014  coes2(7,:)=(/0.0,0.0,0.0,0.0,1.551e-5,-0.002103,0.03268,-0.4304,0.3751/)
1015 
1016  c_ice=(aice(ia,1)+aice(ib,1))/2
1017  floe_diam=200
1018  dis=((vy(ia)-vy(ib))**2+(vx(ia)-vx(ib))**2)**0.5*2.0
1019  if(c_ice>0.00001)then
1020  mean_diam_u=0.0
1021  mean_diam_d=0.0
1022  mean_diam=0.0
1023  DO mm=1,3
1024  mean_diam_u=mean_diam_u+(ee**2*ff)**mm*ee**(-mm)*(floe_diam+floe_diam)/2
1025  mean_diam_d=mean_diam_d+(ee**2*ff)**mm
1026  END DO
1027  mean_diam=mean_diam_u/mean_diam_d
1028  if(c_ice.eq.1)then
1029  icethi=3.0
1030  else
1031  icethi=0.3*(1.0+4*c_ice)
1032  endif
1033  DO iss=1,msc
1034  spcsigl=spcsig(iss)
1035  deploc=dep3(i)
1036 ! CALL KSCIP1(1,SPCSIGL,DEPLOC,KWAVEL,CGOL,NN,ND)
1037  twave=1/(spcsigl/pi2_w)
1038  IF(twave<=16.0)THEN
1039  IF(twave<=6.0)twave=6.0
1040  iceatts(1)=coes2(1,1)*twave**8+coes2(1,2)*twave**7+coes2(1,3)*twave**6+coes2(1,4)*twave**5+coes2(1,5)*twave**4+&
1041  coes2(1,6)*twave**3+coes2(1,7)*twave**2+coes2(1,8)*twave+coes2(1,9)
1042 
1043  iceatts(2)=coes2(2,1)*twave**8+coes2(2,2)*twave**7+coes2(2,3)*twave**6+coes2(2,4)*twave**5+coes2(2,5)*twave**4+&
1044  coes2(2,6)*twave**3+coes2(2,7)*twave**2+coes2(2,8)*twave+coes2(2,9)
1045 
1046  iceatts(3)=coes2(3,1)*twave**8+coes2(3,2)*twave**7+coes2(3,3)*twave**6+coes2(3,4)*twave**5+coes2(3,5)*twave**4+&
1047  coes2(3,6)*twave**3+coes2(3,7)*twave**2+coes2(3,8)*twave+coes2(3,9)
1048 
1049  iceatts(4)=coes2(4,1)*twave**8+coes2(4,2)*twave**7+coes2(4,3)*twave**6+coes2(4,4)*twave**5+coes2(4,5)*twave**4+&
1050  coes2(4,6)*twave**3+coes2(4,7)*twave**2+coes2(4,8)*twave+coes2(4,9)
1051 
1052  iceatts(5)=coes2(5,1)*twave**8+coes2(5,2)*twave**7+coes2(5,3)*twave**6+coes2(5,4)*twave**5+coes2(5,5)*twave**4+&
1053  coes2(5,6)*twave**3+coes2(5,7)*twave**2+coes2(5,8)*twave+coes2(5,9)
1054 
1055  iceatts(6)=coes2(6,1)*twave**8+coes2(6,2)*twave**7+coes2(6,3)*twave**6+coes2(6,4)*twave**5+coes2(6,5)*twave**4+&
1056  coes2(6,6)*twave**3+coes2(6,7)*twave**2+coes2(6,8)*twave+coes2(6,9)
1057 
1058  iceatts(7)=coes2(7,1)*twave**8+coes2(7,2)*twave**7+coes2(7,3)*twave**6+coes2(7,4)*twave**5+coes2(7,5)*twave**4+&
1059  coes2(7,6)*twave**3+coes2(7,7)*twave**2+coes2(7,8)*twave+coes2(7,9)
1060  IF (icethi<=0.4)THEN
1061  iceatt=iceatts(1)
1062  ELSE IF(icethi==0.6)THEN
1063  iceatt=iceatts(2)
1064  ELSE IF(icethi==0.8)THEN
1065  iceatt=iceatts(3)
1066  ELSE IF(icethi==1.2)THEN
1067  iceatt=iceatts(4)
1068  ELSE IF(icethi==1.6)THEN
1069  iceatt=iceatts(5)
1070  ELSE IF(icethi==2.4)THEN
1071  iceatt=iceatts(6)
1072  ELSE IF(icethi>=3.2)THEN
1073  iceatt=iceatts(7)
1074  ELSE IF(icethi>0.4.AND.icethi<0.6)THEN
1075  iceatt1=iceatts(1)
1076  iceatt2=iceatts(2)
1077  rate=(icethi-0.4)/0.2
1078  iceatt=iceatt1*(1-rate)+iceatt2*rate
1079  ELSE IF(icethi>0.6.AND.icethi<0.8)THEN
1080  iceatt1=iceatts(2)
1081  iceatt2=iceatts(3)
1082  rate=(icethi-0.6)/0.2
1083  iceatt=iceatt1*(1-rate)+iceatt2*rate
1084  ELSEIF(icethi>0.8.AND.icethi<1.2)THEN
1085  iceatt1=iceatts(3)
1086  iceatt2=iceatts(4)
1087  rate=(icethi-0.8)/0.4
1088  iceatt=iceatt1*(1-rate)+iceatt2*rate
1089  ELSE IF(icethi>1.2.AND.icethi<1.6)THEN
1090  iceatt1=iceatts(4)
1091  iceatt2=iceatts(5)
1092  rate=(icethi-1.2)/0.4
1093  iceatt=iceatt1*(1-rate)+iceatt2*rate
1094  ELSE IF(icethi>1.6.AND.icethi<2.4)THEN
1095  iceatt1=iceatts(5)
1096  iceatt2=iceatts(6)
1097  rate=(icethi-1.6)/0.8
1098  iceatt=iceatt1*(1-rate)+iceatt2*rate
1099  ELSE IF(icethi>2.4.AND.icethi<3.2)THEN
1100  iceatt1=iceatts(6)
1101  iceatt2=iceatts(7)
1102  rate=(icethi-2.4)/0.8
1103  iceatt=iceatt1*(1-rate)+iceatt2*rate
1104  ENDIF
1105  iceatt=exp(iceatt)
1106  iceatt=iceatt*c_ice/mean_diam
1107  ratio(iss)=exp(-iceatt*dis)
1108  ELSE
1109  iceatt=0.02*exp(-0.386*twave)*c_ice
1110  ratio(iss)=exp(-iceatt*dis)
1111  ENDIF
1112  END DO
1113  else
1114  DO iss=1,msc
1115  ratio(iss)=1.0
1116  end do
1117  endif
1118 
1119  RETURN
1120  END SUBROUTINE actt2
1121 !==============================================================================|
1122 
1123 END MODULE mod_action_ex
vars_wave
Definition: mod_main_wave.f90:14
m_genarr::spcsig
real, dimension(:), allocatable, save spcsig
Definition: swmod2.f90:767
swcomm3::pdiffr
real, dimension(mdiffr) pdiffr
Definition: swmod1.f90:2142
mod_action_ex::action_deallo
subroutine action_deallo
Definition: mod_action_ex.f90:44
swcomm3::icur
integer icur
Definition: swmod1.f90:2125
swcomm3::depmin
real depmin
Definition: swmod1.f90:2133
mod_action_ex::n32
real, dimension(:,:,:), allocatable n32
Definition: mod_action_ex.f90:18
m_genarr::spcdir
real, dimension(:,:), allocatable, save spcdir
Definition: swmod2.f90:767
mod_northpole::node_northpole
integer node_northpole
Definition: mod_northpole.f90:48
ice_state
Definition: ice_state.f90:64
mod_action_ex::adv_hl
subroutine adv_hl(FF1, I, I1, I2)
Definition: mod_action_ex.f90:661
vars_wave::ac2
real, dimension(:,:,:), allocatable ac2
Definition: mod_main_wave.f90:24
m_genarr
Definition: swmod2.f90:665
swcomm3::mdc
integer mdc
Definition: swmod1.f90:1672
swcomm3::pi2_w
real pi2_w
Definition: swmod1.f90:1722
mod_action_ex::adv_hl_fij
subroutine adv_hl_fij(N32_TMPP3, N32_TMPP2, I, IA, IB)
Definition: mod_action_ex.f90:787
mod_northpole
Definition: mod_northpole.f90:40
mod_action_ex::algorithm_fct
subroutine algorithm_fct(CAS, IG, DTW, IDCMIN, IDCMAX)
Definition: mod_action_ex.f90:58
swcomm3::jvx2
integer jvx2
Definition: swmod1.f90:1609
all_vars::vx
real(sp), dimension(:), allocatable, target vx
Definition: mod_main.f90:1001
all_vars::vy
real(sp), dimension(:), allocatable, target vy
Definition: mod_main.f90:1002
mod_action_ex
Definition: mod_action_ex.f90:12
kscip1
subroutine kscip1(MMT, SIG, D, K, CG, N, ND)
Definition: swanser.f90:504
ice_state::aice
real(kind=dbl_kind), dimension(:,:), allocatable, target, save aice
Definition: ice_state.f90:82
swcomm3::jdp2
integer jdp2
Definition: swmod1.f90:1601
mod_action_ex::actt2
subroutine actt2(IA, IB, RATIO)
Definition: mod_action_ex.f90:979
all_vars
Definition: mod_main.f90:959
mod_action_ex::adv_n
subroutine adv_n(DTW)
Definition: mod_action_ex.f90:327
sproxy2
subroutine sproxy2(CAXL, CAYL, CG0L, ECOSL, ESINL, UX2L, UY2L)
Definition: swancom5.f90:260
swcomm3::msc
integer msc
Definition: swmod1.f90:1673
mod_action_ex::tridag
subroutine tridag(A, B, C, R, U, N)
Definition: mod_action_ex.f90:301
kscip2
subroutine kscip2(S, D, CG)
Definition: swanser.f90:650
swcomm3::idiffr
integer idiffr
Definition: swmod1.f90:2132
swcomm3::jvy2
integer jvy2
Definition: swmod1.f90:1610
m_parall
Definition: swmod2.f90:797
mod_action_ex::action_allo
subroutine action_allo
Definition: mod_action_ex.f90:27
mod_action_ex::algorithm_crank_nicolson
subroutine algorithm_crank_nicolson(CAD, IG, DTW, IDCMIN, IDCMAX, DD)
Definition: mod_action_ex.f90:239
mod_action_ex::n31
real, dimension(:,:), allocatable n31
Definition: mod_action_ex.f90:17
vars_wave::compda
real, dimension(:,:), allocatable compda
Definition: mod_main_wave.f90:24
Generated by   doxygen 1.8.15