module mod_coriolis use mod_parameters use mod_variables implicit none contains subroutine convert2jmax implicit none if(SWITCH_coriolis == 1) then jstart=1; jfinal=jmax end if status = 1 end subroutine convert2jmax subroutine coriolis2 implicit none integer :: j, m, il, ir, jm1, jm0, jp1 real(8) :: jj, mm, ff complex(8) :: fcjm1_vjm1, fcjm1_vj complex(8) :: fcj_vjm1, fcj_vj, fcj_vjp1 complex(8) :: fcjp1_vj, fcjp1_vjp1 do il = 1, numoflayers layer(il)%fc(:,:) = dcmplx(0.0d0, 0.0d0) end do do j = 1, jmax jj = dble(j) do m = 0, min(2, j), 2 mm = dble(m) jm1=jml2(j,m, j-1); jm0=jml2(j,m,j); jp1=jml2(j,m,j+1) fcjm1_vjm1 = (Im/jj) * dsqrt( (jj-1.0d0) * (jj**2 - mm**2) / (2.0d0*jj-1.0d0) ) fcjm1_vj = -(Im*mm/jj) fcj_vjm1 = (Im/jj) * dsqrt( (jj+1.0d0) * (jj**2 - mm**2) / (2.0d0*jj+1.0d0) ) fcj_vj = -(Im*mm/(jj*(jj+1.0d0))) fcj_vjp1 = (Im/(jj+1.0d0)) * dsqrt( (jj*((jj+1.0d0)**2 - mm**2)) / (2.0d0*jj+1.0d0) ) fcjp1_vj = (Im*mm)/(jj+1.0d0) fcjp1_vjp1 = (Im/(jj+1.0d0)) * dsqrt( ((jj+2.0d0)*((jj+1.0d0)**2 - mm**2)) / (2.0d0*jj+3.0d0) ) do il = layer_ocean, layer_ocean do ir = 1, layer(il)%n-1 ff = -2.0d0 * layer(il)%rho_l(ir+1) * ang_vel layer(il)%fc(ir, jm1) = ff*(fcjm1_vjm1 * layer(il)%vel(ir+1, jml2(j-1,m, j-1)) + fcjm1_vj * layer(il)%vel(ir+1, jml2(j ,m,j-1))) layer(il)%fc(ir, jm0) = ff*(fcj_vjm1 * layer(il)%vel(ir+1, jml2(j-1,m, j )) + fcj_vj * layer(il)%vel(ir+1, jml2(j ,m,j )) + fcj_vjp1 * layer(il)%vel(ir+1, jml2(j+1, m, j))) layer(il)%fc(ir, jp1) = ff*(fcjp1_vj * layer(il)%vel(ir+1, jml2(j ,m, j+1)) + fcjp1_vjp1 * layer(il)%vel(ir+1, jml2(j+1 ,m,j+1))) end do end do end do end do end subroutine coriolis2 subroutine coriolis implicit none real(8) :: cc, c6, co, coef1, fco integer :: j2, deg, ord, l, jml, jml1, i, k do i = 1, numoflayers layer(i)%fc1(:,:) = dcmplx(0.0d0, 0.0d0) end do do deg = 0, jmax order: do ord = 0, min(2,deg), 2 vec_l: do l = abs(deg - 1), deg + 1 jml = jml2(deg, ord, l) coef1 = (-1.0d0)**(deg + l) ddeg_2: do j2 = max(deg-1, 1, ord, l-1), min(deg+1, jmax, l+1) jml1 = jml2(j2, ord, l) call cleb_1(j2, ord, 1, 0, deg, ord, cc) call six_j(j2, 1, l, deg, 1, c6) do i = layer_ocean, layer_ocean do k = 1, layer(i)%n-1 co = 2.0d0 * dsqrt(6.0d0) * ang_vel * layer(i)%rho_l(k+1) fco = co * coef1 * dsqrt(dble(2*j2 + 1)) * cc * c6 layer(i)%fc1(k, jml) = layer(i)%fc1(k, jml) + (Im * fco * layer(i)%vel(k+1, jml1)) end do end do end do ddeg_2 end do vec_l end do order end do end subroutine coriolis subroutine coriolis_pvel implicit none real(8) :: cc, c6, co, coef1, fco integer :: j2, deg, ord, l, jml, jml1, i, k do i = 1, numoflayers layer(i)%fc0(:,:) = dcmplx(0.0d0, 0.0d0) end do do deg = 0, jmax order: do ord = 0, min(2,deg), 2 vec_l: do l = abs(deg - 1), deg + 1 jml = jml2(deg, ord, l) coef1 = (-1.0d0)**(deg + l) ddeg_2: do j2 = max(deg-1, 1, ord, l-1), min(deg+1, jmax, l+1) jml1 = jml2(j2, ord, l) call cleb_1(j2, ord, 1, 0, deg, ord, cc) call six_j(j2, 1, l, deg, 1, c6) do i = layer_ocean, layer_ocean do k = 1, layer(i)%n-1 co = 2.0d0 * dsqrt(6.0d0) * ang_vel * layer(i)%rho_l(k+1) fco = co * coef1 * dsqrt(dble(2*j2 + 1)) * cc * c6 layer(i)%fc0(k, jml) = layer(i)%fc0(k, jml) + (Im * fco * layer(i)%pvel(k+1, jml1)) end do end do end do ddeg_2 end do vec_l end do order end do end subroutine coriolis_pvel subroutine six_j(j1, j2, j, l2, l, c6) ! 6-j symbol : {j1,j2,j ! 1,l2,l} !..................................................................... implicit none integer :: j1, j2, j, l2, l real(8) :: c6, a, b, c, s, zn real(8) :: hh, dd c6=0.0d0 if(l2.gt.(j+1).or.iabs(j-1).gt.l2) return if(j.gt.(j1+j2).or.iabs(j1-j2).gt.j) return if(l.gt.(j2+1).or.iabs(j2-1).gt.l) return a=dble(j1) b=dble(j2) c=dble(j) s=a+b+c zn=(-1.0d0)**(j1+j2+j) if((j-l2).lt.0) then if((j2-l).lt.0)then hh=(s+2)*(s+3)*(s-a-a+1)*(s-a-a+2) dd=(b+b+1)*(b+1)*(b+b+3)*(c+c+1)*(c+1)*(c+c+3) end if if((j2-l).eq.0)then hh=(s+2)*(s-c-c)*(s-b-b+1)*(s-a-a+1) dd=b*(b+b+1)*(b+1)*(c+c+1)*(c+1)*(c+c+3) zn=-zn end if if((j2-l).gt.0)then hh=(s-c-c-1)*(s-c-c)*(s-b-b+1)*(s-b-b+2) dd=(b+b-1)*b*(b+b+1)*(c+c+1)*(c+1)*(c+c+3) end if end if if((j-l2).eq.0) then if((j2-l).lt.0)then hh=(s+2)*(s-c-c+1)*(s-b-b)*(s-a-a+1) dd=(b+b+1)*(b+1)*(b+b+3)*c*(c+c+1)*(c+1) zn=-zn end if if((j2-l).eq.0)then hh=-a*(a+1)+b*(b+1)+c*(c+1) dd=b*(b+b+1)*(b+1)*c*(c+c+1)*(c+1) c6=-zn*hh/dsqrt(dd)/2.d00; return end if if((j2-l).gt.0)then hh=(s+1)*(s-c-c)*(s-b-b+1)*(s-a-a) dd=(b+b-1)*b*(b+b+1)*c*(c+c+1)*(c+1) end if end if if((j-l2).gt.0) then if((j2-l).lt.0)then hh=(s-c-c+1)*(s-c-c+2)*(s-b-b-1)*(s-b-b) dd=(b+b+1)*(b+1)*(b+b+3)*(c+c-1)*c*(c+c+1) end if if((j2-l).eq.0)then hh=(s+1)*(s-c-c+1)*(s-b-b)*(s-a-a) dd=b*(b+b+1)*(b+1)*(c+c-1)*c*(c+c+1) end if if((j2-l).gt.0)then hh=s*(s+1)*(s-a-a-1)*(s-a-a) dd=(b+b-1)*b*(b+b+1)*(c+c-1)*c*(c+c+1) end if end if c6=zn*dsqrt(hh/dd)*0.5d0 end subroutine six_j subroutine cleb_1(j1,m1,j2,m2,j,m,cg) ! ! j m ! C ! j1 m1 1 m2 !........................................................................ implicit none integer :: j1, m1, j2, m2, j, m real(8) :: cg, c, g, zn cg=0.0d0 if(j2.ne.1) return if(abs(j1-j).gt.1.or.(j1+j).eq.0) return if(abs(m2).gt.1.or.abs(m1).gt.j1) return if(m.ne.(m1+m2)) return c=dble(j) g=dble(m) zn=1.0d0 if(m2.lt.0)then if((j1-j).lt.0)then cg=(c-g-1.0d0)*(c-g)/(c+c-1.0d0)/(c+c) elseif((j1-j).eq.0)then cg=(c+g+1.0d0)*(c-g)/(c+1.0d0)/(c+c) else cg=(c+g+2.0d0)*(c+g+1.0d0)/(c+c+2.0d0)/(c+c+3.0d0) end if end if if(m2.eq.0)then if((j1-j).lt.0)then cg=(c+g)*(c-g)/(c+c-1.0d0)/c elseif((j1-j).eq.0)then cg=g/dsqrt(c*(c+1.0d0)); return else cg=(c+g+1.0d0)*(c-g+1.0d0)/(c+1.0d0)/(c+c+3.0d0) zn=-1.0d0 end if end if if(m2.gt.0)then if((j1-j).lt.0)then cg=(c+g-1.0d0)*(c+g)/(c+c-1.0d0)/(c+c) elseif((j1-j).eq.0)then cg=(c+g)*(c-g+1.0d0)/(c+1.0d0)/(c+c) zn=-1.0d0 else cg=(c-g+1.0d0)*(c-g+2.0d0)/(c+c+2.0d0)/(c+c+3.0d0) end if end if cg=zn*dsqrt(cg) end subroutine cleb_1 end module mod_coriolis