mod_spherical.f90

Go to the documentation of this file.

!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

!==============================================================================|
!  MODULE CONTAINING SUBROUTINES USED TO SET UP SPHERICAL COORDINATE SYSTEM    |
!  FLUXES                                                                      |
!==============================================================================|

MODULE mod_spherical
  USE control
  USE mod_prec
  IMPLICIT NONE
  SAVE
  REAL(sp), ALLOCATABLE :: dltxne(:,:)
  REAL(sp), ALLOCATABLE :: dltyne(:,:)

  REAL(sp), ALLOCATABLE :: deltux(:,:)
  REAL(sp), ALLOCATABLE :: deltuy(:,:)
  REAL(sp), ALLOCATABLE :: sitau(:,:)


  INTERFACE arc
     MODULE PROCEDURE arc_flt
     MODULE PROCEDURE arc_dbl
  END INTERFACE

  INTERFACE arcc
     MODULE PROCEDURE arcc_flt
     MODULE PROCEDURE arcc_dbl
  END INTERFACE

  INTERFACE area
     MODULE PROCEDURE area_flt
     MODULE PROCEDURE area_dbl
  END INTERFACE

  INTERFACE arcx
     MODULE PROCEDURE arcx_flt
     MODULE PROCEDURE arcx_dbl
  END INTERFACE

  INTERFACE arcy
     MODULE PROCEDURE arcy_flt
     MODULE PROCEDURE arcy_dbl
  END INTERFACE

  INTERFACE arcx_back
     MODULE PROCEDURE arcx_back_flt
     MODULE PROCEDURE arcx_back_dbl
  END INTERFACE


  !===================================================================================|
CONTAINS   !!INCLUDED SUBROUTINES FOLLOW
  !===================================================================================|

  SUBROUTINE arc_dbl(XX1,YY1,XX2,YY2,ARCL)
    !----------------------------------------------------------------------------
    !      function:
    !           calculate the arc lenth for given two point on the spherical plane
    !      input:
    !           xx1,yy1,xx2,yy2 :are longitude and latitude of two points
    !      output:
    !           arcl :  arc lenth of two points in spherical plane
    !-----------------------------------------------------------------------------       

    !  solve the arc length through the earth center
    IMPLICIT NONE
    REAL(DP) :: X1,Y1,X2,Y2,XA,YA,ZA,XB,YB,ZB,AB,AOB
    REAL(DP),INTENT(OUT) :: ARCL
    REAL(DP),INTENT(IN) :: XX1,YY1,XX2,YY2

    x1=xx1*deg2rad
    y1=yy1*deg2rad

    x2=xx2*deg2rad
    y2=yy2*deg2rad

    ! USE DOUBLE PRECISION COS AND SIN
    xa=dcos(y1)*dcos(x1)
    ya=dcos(y1)*dsin(x1)
    za=dsin(y1)

    xb=dcos(y2)*dcos(x2)
    yb=dcos(y2)*dsin(x2)
    zb=dsin(y2)

    ab=dsqrt((xb-xa)**2+(yb-ya)**2+(zb-za)**2)
    aob=(2.0_dp -ab*ab)/2.0_dp
    aob=dacos(aob)
    arcl=rearth*aob

    RETURN
  END SUBROUTINE arc_dbl

  SUBROUTINE arc_flt(XX1,YY1,XX2,YY2,ARCL)
    IMPLICIT NONE
    REAL(SPA), INTENT(IN) :: XX1,YY1,XX2,YY2
    REAL(SPA), INTENT(OUT) :: ARCL
    REAL(DP) ARCL_DP
    CALL arc_dbl(dble(xx1),dble(yy1),dble(xx2),dble(yy2),arcl_dp)
    arcl = arcl_dp
  END SUBROUTINE arc_flt


  SUBROUTINE area_dbl(SIDEA,SIDEB,SIDEC,AREA1)
    !--------------------------------------------------------------------
    !      function:
    !           calculate the area of a triangle on a spherical plane
    !      input:
    !           side1,side2 and side3: are 3 arc lenth for one triangle
    !      output:
    !           areal: is area of a triangle on a spherical plane
    !--------------------------------------------------------------------
    IMPLICIT NONE
    REAL(DP), INTENT(IN) :: SIDEA,SIDEB,SIDEC
    REAL(DP), INTENT(OUT) :: AREA1
    REAL(DP) :: SIDE1,SIDE2,SIDE3
    REAL(DP) :: PSUM,PM,QMJC

    side1=sidea/rearth
    side2=sideb/rearth
    side3=sidec/rearth

    ! SLOWER TO CHECK THEN TO CALCULATE
    !   IF(SIDE1 == 0.0_DP .OR. SIDE2 == 0.0_DP .OR. SIDE3 == 0.0_DP)THEN
    !     AREA1=0.0_DP
    !   ELSE

    psum=0.5_dp*(side1+side2+side3)
    pm=dsin(psum)*dsin(psum-side1)*dsin(psum-side2)*dsin(psum-side3)
    pm=dsqrt(pm)/(2.0_dp*dcos(side1*0.5_dp)*dcos(side2*0.5_dp)*dcos(side3*0.5_dp))
    qmjc = 2.0_dp*dasin(pm)

    area1=rearth*rearth*qmjc

    !   END IF

    RETURN
  END SUBROUTINE area_dbl

  SUBROUTINE area_flt(SIDE1,SIDE2,SIDE3,AREA1)
    IMPLICIT NONE
    REAL(SPA), INTENT(IN) :: SIDE1,SIDE2,SIDE3
    REAL(SPA), INTENT(OUT) :: AREA1
    REAL(DP) :: AREA_DP

    CALL area_dbl(dble(side1),dble(side2),dble(side3),area_dp)
    area1=area_dp

  END SUBROUTINE area_flt


  SUBROUTINE arcc_dbl(XX1,YY1,XX2,YY2,XXC,YYC)
    IMPLICIT NONE
    REAL(DP), INTENT(OUT) :: XXC,YYC
    REAL(DP), INTENT(IN) :: XX1,YY1,XX2,YY2
    REAL(DP) :: X1,Y1,X2,Y2

    x1=xx1*deg2rad
    y1=yy1*deg2rad

    x2=xx2*deg2rad
    y2=yy2*deg2rad

    xxc=dcos(y1)*dsin(x1)+dcos(y2)*dsin(x2)
    !   XXC=XXC/(COS(Y1)*COS(X1)+COS(Y2)*COS(X2))
    !   XXC=ATAN(XXC)
    xxc=datan2(xxc,(dcos(y1)*dcos(x1)+dcos(y2)*dcos(x2)))
    xxc=xxc/deg2rad

    !   IF(XXC .LT. 0.0) XXC=180.0+XXC
    IF(xxc < 0.0_dp) xxc=360.0_dp+xxc

    yyc=dcos(y1)*dcos(y1)+dcos(y2)*dcos(y2)+2.0_dp*dcos(y1)*dcos(y2)*dcos(x1-x2)
    !   YYC=SQRT(YYC)/(SIN(Y1)+SIN(Y2))
    yyc=datan2(dsqrt(yyc),(dsin(y1)+dsin(y2)))
    !   YYC=ATAN(YYC)
    yyc=90.0_dp-yyc/deg2rad

    RETURN
  END SUBROUTINE arcc_dbl

  SUBROUTINE arcc_flt(XX1,YY1,XX2,YY2,XXC,YYC)
    IMPLICIT NONE
    REAL(SPA) :: XX1,YY1,XX2,YY2
    REAL(SPA) :: XXC,YYC
    REAL(DP) :: XXC_DP,YYC_DP

    CALL arcc_dbl(dble(xx1),dble(yy1),dble(xx2),dble(yy2),xxc_dp,yyc_dp)
    xxc = xxc_dp
    yyc = yyc_dp

  END SUBROUTINE arcc_flt


  SUBROUTINE arcx_dbl(XX1,YY1,XX2,YY2,ARCX1)
    IMPLICIT NONE
    REAL(DP), INTENT(IN) :: XX1,YY1,XX2,YY2
    REAL(DP), INTENT(OUT)::ARCX1

    REAL(DP) :: X1,Y1,X2,Y2,TY
    REAL(DP) :: XTMP          

    IF(xx1 == xx2)THEN
       arcx1=0.0_dp
    ELSE
       x1=xx1*deg2rad
       y1=yy1*deg2rad

       x2=xx2*deg2rad
       y2=yy2*deg2rad

       xtmp  = x2-x1
       IF(xtmp >  pi)THEN
          xtmp = real(-2*pi,dp)+xtmp
       ELSE IF(xtmp < -pi)THEN
          xtmp =  real(2*pi,dp)+xtmp
       END IF

       ty=0.5_dp*(y2+y1)
       arcx1=rearth*dcos(ty)*xtmp
    END IF

    RETURN
  END SUBROUTINE arcx_dbl

  SUBROUTINE arcy_dbl(XX1,YY1,XX2,YY2,ARCY1)
    IMPLICIT NONE
    REAL(DP), INTENT(IN) :: XX1,YY1,XX2,YY2
    REAL(DP), INTENT(OUT)::ARCY1

    REAL(DP) :: X1,Y1,X2,Y2,TY
    REAL(DP) :: YTMP          

    IF(yy1 == yy2)THEN
       arcy1=0.0_dp
    ELSE
       x1=xx1*deg2rad
       y1=yy1*deg2rad

       x2=xx2*deg2rad
       y2=yy2*deg2rad

       ytmp  = y2-y1
       IF(ytmp >  pi)THEN
          ytmp = real(-2*pi,dp)+ytmp
       ELSE IF(ytmp < -pi)THEN
          ytmp =  real(2*pi,dp)+ytmp
       END IF

       arcy1=rearth*ytmp
    END IF

    RETURN
  END SUBROUTINE arcy_dbl

  SUBROUTINE arcy_flt(XX1,YY1,XX2,YY2,ARCY1)
    IMPLICIT NONE
    REAL(SPA), INTENT(IN) :: XX1,YY1,XX2,YY2
    REAL(SPA), INTENT(OUT)::ARCY1
    
    REAL(DP) ::ARCY_DP

    CALL arcy_dbl(dble(xx1),dble(yy1),dble(xx2),dble(yy2),arcy_dp)
    arcy1 = arcy_dp

  END SUBROUTINE arcy_flt

  SUBROUTINE arcx_flt(XX1,YY1,XX2,YY2,ARCX1)
    IMPLICIT NONE
    REAL(SPA), INTENT(IN) :: XX1,YY1,XX2,YY2
    REAL(SPA), INTENT(OUT)::ARCX1
    
    REAL(DP) ::ARCX_DP

    CALL arcx_dbl(dble(xx1),dble(yy1),dble(xx2),dble(yy2),arcx_dp)
    arcx1 = arcx_dp

  END SUBROUTINE arcx_flt

  SUBROUTINE arcx_back_dbl(XX1,YY1,XX2,YY2,ARCX1)
    IMPLICIT NONE
    REAL(DP), INTENT(IN) :: XX1,YY1,XX2,YY2
    REAL(DP), INTENT(OUT) :: ARCX1

    INTEGER I
    INTEGER,PARAMETER ::NX=500
    REAL(DP) :: X1,Y1,X2,Y2,TY,A1,A2,B1,B2,C1,C2,A,B,C,X(NX+1),Y(NX+1)
    REAL(DP) :: XTMP          

    IF(xx1 == xx2)THEN
       arcx1=0.
    ELSE
       x1=xx1*deg2rad
       y1=yy1*deg2rad

       x2=xx2*deg2rad
       y2=yy2*deg2rad

       x(1)=x1
       y(1)=y1
       x(nx+1)=x2
       y(nx+1)=y2

       xtmp=x(nx+1)-x(1)
       IF(xtmp >  pi)THEN
          xtmp = real(-2*pi,dp)+xtmp
       ELSE IF(xtmp < -pi)THEN
          xtmp =  real(2*pi,dp)+xtmp
       END IF

       DO i=2,nx
          x(i)=x(i-1)+xtmp/float(nx)
          !       x(i)=x(i-1)+(x(nx+1)-x(1))/float(nx)
       END DO

       a1=dcos(y(1))*dcos(x(1))
       a2=dcos(y(nx+1))*dcos(x(nx+1))

       b1=dcos(y(1))*dsin(x(1))
       b2=dcos(y(nx+1))*dsin(x(nx+1))

       c1=dsin(y(1))
       c2=dsin(y(nx+1))

       a=a1*b2-a2*b1
       b=b1*c2-b2*c1
       c=a2*c1-a1*c2

       DO i=2,nx
          y(i)=-b*dcos(x(i))-c*dsin(x(i))
          y(i)=y(i)/a
          y(i)=datan(y(i))
       END DO

       arcx1=0.
       DO i=1,nx
          ty=0.5*(y(i)+y(i+1))
          xtmp=x(i+1)-x(i)
          IF(xtmp >  pi)THEN
             xtmp = real(-2*pi,dp)+xtmp
          ELSE IF(xtmp < -pi)THEN
             xtmp =  real(2*pi,dp)+xtmp
          END IF
          arcx1=arcx1+rearth*dcos(ty)*xtmp
          !       arcx1=arcx1+rearth*cos(ty)*(x(i+1)-x(i))
       END DO
    END IF

    RETURN
  END SUBROUTINE arcx_back_dbl

    SUBROUTINE arcx_back_flt(XX1,YY1,XX2,YY2,ARCX1)
    IMPLICIT NONE
    REAL(SPA), INTENT(IN) :: XX1,YY1,XX2,YY2
    REAL(SPA), INTENT(OUT)::ARCX1
    
    REAL(DP) ::ARCX_DP

    CALL arcx_back_dbl(dble(xx1),dble(yy1),dble(xx2),dble(yy2),arcx_dp)
    arcx1 = arcx_dp

  END SUBROUTINE arcx_back_flt

  SUBROUTINE alloc_sphere_vars
    USE lims
    INTEGER NCT
    INTEGER NDB
    ndb = 1       !!GWC BASE THIS ON KIND

    nct = nt*3
    ALLOCATE(dltxne(nct,2))       ;dltxne    = zero
    ALLOCATE(dltyne(nct,2))       ;dltyne    = zero
    ALLOCATE(deltux(nt,3))        ;deltux    = zero
    ALLOCATE(deltuy(nt,3))        ;deltuy    = zero
    ALLOCATE(sitau(nt,3))         ;sitau     = zero

    memcnt = memcnt + nct*4*ndb + nt*9*ndb
    RETURN
  END SUBROUTINE alloc_sphere_vars


END MODULE  mod_spherical