Skip to content

Commit

Permalink
CUDA Fortran SCC-HOIST: Loki generated and copied to ecwam/src
Browse files Browse the repository at this point in the history
  • Loading branch information
@MichaelSt98
MichaelSt98 committed Mar 21, 2024
1 parent 00c73c4 commit 0f522dd
Show file tree
Hide file tree
Showing 49 changed files with 9,999 additions and 0 deletions.
182 changes: 182 additions & 0 deletions src/phys-scc-cuf-hoist/airsea.cuf_hoist_new.F90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,182 @@
! (C) Copyright 1989- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!
MODULE AIRSEA_CUF_HOIST_NEW_MOD
CONTAINS
ATTRIBUTES(DEVICE) SUBROUTINE AIRSEA_CUF_HOIST_NEW (KIJS, KIJL, HALP, U10, U10DIR, TAUW, TAUWDIR, RNFAC, US, Z0, Z0B, CHRNCK, &
& ICODE_WND, IUSFG, ACD, ALPHA, ALPHAMAX, ALPHAMIN, ANG_GC_A, ANG_GC_B, ANG_GC_C, BCD, BETAMAXOXKAPPA2, BMAXOKAP, &
& C2OSQRTVG_GC, CDMAX, CHNKMIN_U, CM_GC, DELKCC_GC_NS, DELKCC_OMXKM3_GC, EPS1, EPSMIN, EPSUS, G, GM1, LLCAPCHNK, LLGCBZ0, &
& LLNORMAGAM, NWAV_GC, OM3GMKM_GC, OMXKM3_GC, RN1_RN, RNU, RNUM, SQRTGOSURFT, WSPMIN, XKAPPA, XKMSQRTVGOC2_GC, XKM_GC, XK_GC, &
& XLOGKRATIOM1_GC, XNLEV, ZALP, ICHNK, NCHNK, IJ)

! ----------------------------------------------------------------------

!**** *AIRSEA* - DETERMINE TOTAL STRESS IN SURFACE LAYER.

! P.A.E.M. JANSSEN KNMI AUGUST 1990
! JEAN BIDLOT ECMWF FEBRUARY 1999 : TAUT is already
! SQRT(TAUT)
! JEAN BIDLOT ECMWF OCTOBER 2004: QUADRATIC STEP FOR
! TAUW

!* PURPOSE.
! --------

! COMPUTE TOTAL STRESS.

!** INTERFACE.
! ----------

! *CALL* *AIRSEA (KIJS, KIJL, FL1, WAVNUM,
! HALP, U10, U10DIR, TAUW, TAUWDIR, RNFAC,
! US, Z0, Z0B, CHRNCK, ICODE_WND, IUSFG)*

! *KIJS* - INDEX OF FIRST GRIDPOINT.
! *KIJL* - INDEX OF LAST GRIDPOINT.
! *FL1* - SPECTRA
! *WAVNUM* - WAVE NUMBER
! *HALP* - 1/2 PHILLIPS PARAMETER
! *U10* - WINDSPEED U10.
! *U10DIR* - WINDSPEED DIRECTION.
! *TAUW* - WAVE STRESS.
! *TAUWDIR* - WAVE STRESS DIRECTION.
! *RNFAC* - WIND DEPENDENT FACTOR USED IN THE GROWTH RENORMALISATION.
! *US* - OUTPUT OR OUTPUT BLOCK OF FRICTION VELOCITY.
! *Z0* - OUTPUT BLOCK OF ROUGHNESS LENGTH.
! *Z0B* - BACKGROUND ROUGHNESS LENGTH.
! *CHRNCK* - CHARNOCK COEFFICIENT
! *ICODE_WND* SPECIFIES WHICH OF U10 OR US HAS BEEN FILED UPDATED:
! U10: ICODE_WND=3 --> US will be updated
! US: ICODE_WND=1 OR 2 --> U10 will be updated
! *IUSFG* - IF = 1 THEN USE THE FRICTION VELOCITY (US) AS FIRST GUESS in TAUT_Z0
! 0 DO NOT USE THE FIELD US


! ----------------------------------------------------------------------

USE Z0WAVE_CUF_HOIST_NEW_MOD, ONLY: Z0WAVE_CUF_HOIST_NEW
USE TAUT_Z0_CUF_HOIST_NEW_MOD, ONLY: TAUT_Z0_CUF_HOIST_NEW
USE cudafor
USE PARKIND_WAVE, ONLY: JWIM, JWRU, JWRB

USE YOWPARAM, ONLY: NFRE, NANG
USE YOWTEST, ONLY: IU06


! ----------------------------------------------------------------------
IMPLICIT NONE
INTEGER(KIND=JWIM), VALUE, INTENT(IN) :: KIJS
INTEGER(KIND=JWIM), VALUE, INTENT(IN) :: KIJL
INTEGER(KIND=JWIM), VALUE, INTENT(IN) :: ICODE_WND
INTEGER(KIND=JWIM), VALUE, INTENT(IN) :: IUSFG
REAL(KIND=JWRB), INTENT(IN) :: HALP(KIJL)
REAL(KIND=JWRB), INTENT(IN) :: RNFAC(KIJL)
REAL(KIND=JWRB), INTENT(IN) :: U10DIR(KIJL, NCHNK)
REAL(KIND=JWRB), INTENT(IN) :: TAUW(KIJL, NCHNK)
REAL(KIND=JWRB), INTENT(IN) :: TAUWDIR(KIJL, NCHNK)
REAL(KIND=JWRB), INTENT(INOUT) :: U10(KIJL, NCHNK)
REAL(KIND=JWRB), INTENT(INOUT) :: US(KIJL, NCHNK)
REAL(KIND=JWRB), INTENT(OUT) :: Z0(KIJL, NCHNK)
REAL(KIND=JWRB), INTENT(OUT) :: Z0B(KIJL, NCHNK)
REAL(KIND=JWRB), INTENT(OUT) :: CHRNCK(KIJL, NCHNK)

INTEGER(KIND=JWIM), VALUE, INTENT(IN) :: IJ
INTEGER(KIND=JWIM) :: I
INTEGER(KIND=JWIM) :: J

REAL(KIND=JWRB) :: XI
REAL(KIND=JWRB) :: XJ
REAL(KIND=JWRB) :: DELI1
REAL(KIND=JWRB) :: DELI2
REAL(KIND=JWRB) :: DELJ1
REAL(KIND=JWRB) :: DELJ2
REAL(KIND=JWRB) :: UST2
REAL(KIND=JWRB) :: ARG
REAL(KIND=JWRB) :: SQRTCDM1
REAL(KIND=JWRB) :: XKAPPAD
REAL(KIND=JWRB) :: XLOGLEV
REAL(KIND=JWRB) :: XLEV
REAL(KIND=JWRB), VALUE, INTENT(IN) :: ACD
REAL(KIND=JWRB), VALUE, INTENT(IN) :: ALPHA
REAL(KIND=JWRB), VALUE, INTENT(IN) :: ALPHAMAX
REAL(KIND=JWRB), VALUE, INTENT(IN) :: ALPHAMIN
REAL(KIND=JWRB), VALUE, INTENT(IN) :: ANG_GC_A
REAL(KIND=JWRB), VALUE, INTENT(IN) :: ANG_GC_B
REAL(KIND=JWRB), VALUE, INTENT(IN) :: ANG_GC_C
REAL(KIND=JWRB), VALUE, INTENT(IN) :: BCD
REAL(KIND=JWRB), VALUE, INTENT(IN) :: BETAMAXOXKAPPA2
REAL(KIND=JWRB), VALUE, INTENT(IN) :: BMAXOKAP
REAL(KIND=JWRB), INTENT(IN), DEVICE :: C2OSQRTVG_GC(NWAV_GC)
REAL(KIND=JWRB), VALUE, INTENT(IN) :: CDMAX
REAL(KIND=JWRB), VALUE, INTENT(IN) :: CHNKMIN_U
REAL(KIND=JWRB), INTENT(IN), DEVICE :: CM_GC(NWAV_GC)
REAL(KIND=JWRB), INTENT(IN), DEVICE :: DELKCC_GC_NS(NWAV_GC)
REAL(KIND=JWRB), INTENT(IN), DEVICE :: DELKCC_OMXKM3_GC(NWAV_GC)
REAL(KIND=JWRB), VALUE, INTENT(IN) :: EPS1
REAL(KIND=JWRB), VALUE, INTENT(IN) :: EPSMIN
REAL(KIND=JWRB), VALUE, INTENT(IN) :: EPSUS
REAL(KIND=JWRB), VALUE, INTENT(IN) :: G
REAL(KIND=JWRB), VALUE, INTENT(IN) :: GM1
LOGICAL, VALUE, INTENT(IN) :: LLCAPCHNK
LOGICAL, VALUE, INTENT(IN) :: LLGCBZ0
LOGICAL, VALUE, INTENT(IN) :: LLNORMAGAM
INTEGER(KIND=JWIM), VALUE, INTENT(IN) :: NWAV_GC
REAL(KIND=JWRB), INTENT(IN), DEVICE :: OM3GMKM_GC(NWAV_GC)
REAL(KIND=JWRB), INTENT(IN), DEVICE :: OMXKM3_GC(NWAV_GC)
REAL(KIND=JWRB), VALUE, INTENT(IN) :: RN1_RN
REAL(KIND=JWRB), VALUE, INTENT(IN) :: RNU
REAL(KIND=JWRB), VALUE, INTENT(IN) :: RNUM
REAL(KIND=JWRB), VALUE, INTENT(IN) :: SQRTGOSURFT
REAL(KIND=JWRB), VALUE, INTENT(IN) :: WSPMIN
REAL(KIND=JWRB), VALUE, INTENT(IN) :: XKAPPA
REAL(KIND=JWRB), INTENT(IN), DEVICE :: XKMSQRTVGOC2_GC(NWAV_GC)
REAL(KIND=JWRB), INTENT(IN), DEVICE :: XKM_GC(NWAV_GC)
REAL(KIND=JWRB), INTENT(IN), DEVICE :: XK_GC(NWAV_GC)
REAL(KIND=JWRB), VALUE, INTENT(IN) :: XLOGKRATIOM1_GC
REAL(KIND=JWRB), VALUE, INTENT(IN) :: XNLEV
REAL(KIND=JWRB), VALUE, INTENT(IN) :: ZALP
INTEGER(KIND=JWIM), VALUE, INTENT(IN) :: ICHNK
INTEGER, VALUE, INTENT(IN) :: NCHNK

! ----------------------------------------------------------------------

!* 2. DETERMINE TOTAL STRESS (if needed)
! ----------------------------------

IF (ICODE_WND == 3) THEN

CALL TAUT_Z0_CUF_HOIST_NEW(KIJS, KIJL, IUSFG, HALP(:), U10(:, :), U10DIR(:, :), TAUW(:, :), TAUWDIR(:, :), RNFAC(:), &
& US(:, :), Z0(:, :), Z0B(:, :), CHRNCK(:, :), ACD, ALPHA, ALPHAMAX, ALPHAMIN, ANG_GC_A, ANG_GC_B, ANG_GC_C, BCD, &
& BETAMAXOXKAPPA2, BMAXOKAP, C2OSQRTVG_GC(:), CDMAX, CHNKMIN_U, CM_GC(:), DELKCC_GC_NS(:), DELKCC_OMXKM3_GC(:), EPS1, &
& EPSMIN, EPSUS, G, GM1, LLCAPCHNK, LLGCBZ0, LLNORMAGAM, NWAV_GC, OM3GMKM_GC(:), OMXKM3_GC(:), RN1_RN, RNU, RNUM, &
& SQRTGOSURFT, XKAPPA, XKMSQRTVGOC2_GC(:), XKM_GC(:), XK_GC(:), XLOGKRATIOM1_GC, XNLEV, ZALP, ICHNK, NCHNK, IJ)

ELSE IF (ICODE_WND == 1 .or. ICODE_WND == 2) THEN

!* 3. DETERMINE ROUGHNESS LENGTH (if needed).
! ---------------------------

CALL Z0WAVE_CUF_HOIST_NEW(KIJS, KIJL, US(:, :), TAUW(:, :), U10(:, :), Z0(:, :), Z0B(:, :), CHRNCK(:, :), ALPHA, ALPHAMIN, &
& CHNKMIN_U, EPS1, G, GM1, LLCAPCHNK, ICHNK, NCHNK, IJ)

!* 3. DETERMINE U10 (if needed).
! ---------------------------

XKAPPAD = 1.0_JWRB / XKAPPA
XLOGLEV = LOG(XNLEV)


U10(IJ, ICHNK) = XKAPPAD*US(IJ, ICHNK)*(XLOGLEV - LOG(Z0(IJ, ICHNK)))
U10(IJ, ICHNK) = MAX(U10(IJ, ICHNK), WSPMIN)


END IF


END SUBROUTINE AIRSEA_CUF_HOIST_NEW
END MODULE AIRSEA_CUF_HOIST_NEW_MOD
118 changes: 118 additions & 0 deletions src/phys-scc-cuf-hoist/aki_ice.cuf_hoist_new.F90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,118 @@
! (C) Copyright 1989- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!
MODULE AKI_ICE_CUF_HOIST_NEW_MOD
CONTAINS
ATTRIBUTES(DEVICE) FUNCTION AKI_ICE_CUF_HOIST_NEW (G, XK, DEPTH, RHOW, CITH)

! ----------------------------------------------------------------------

!**** *AKI_ICE* - FUNCTION TO COMPUTE WAVE NUMBER UNDER THE ICE.
! FOR A GIVEN OPEN WATER WAVE NUMBER.

!* PURPOSE.
! -------

! *AKI_ICE* COMPUTES THE REAL WAVE NUMBER UNDER THE ICE AS FUNCTION OF
! OPEN OCEAN WAVE NUMBER, THE WATER DEPTH AND THE ICE THICKNESS.

!** INTERFACE.
! ----------

! *FUNCTION* *AKI_ICE (G,XK,DEPTH,RHOW,CITH))*
! *G* - ACCELERATION OF GRAVITY (m/s**2).
! *XK* - OPEN OCEAN WAVE NUMBER (1/m).
! *DEPTH* - WATER DEPTH (m).
! *RHOW* - WATER DENSITY (kg/m**3).
! *CITH* - SEA ICE THICKNESS (m).

! METHOD.
! -------

! NEWTONS METHOD TO SOLVE THE LINEAR DISPERSION RELATION IN
! SHALLOW WATER UNDER AN INFINITELY LONG ELASTIC FLOATING PLATE
! REPRESENTING THE SEA ICE. THE MECHANICAL PROPERTIES OF THE SEA
! ICE IS GIVEn BY ITS YOUNG MODULUS, THE POISSON'S RATIO AND ITS
! DENSITY (these are fixed for now, see below).

! IF F(x)=0, then solve iteratively
! x(n+1) = x(n) - F(x(n))/F'(x(n))

! WHERE F'(x) IS THE FIRST DERIVATIVE OF F WITH RESPECT TO x.

! EXTERNALS.
! ----------

! NONE.

! REFERENCE.
! ----------

! FOX AND SQUIRE, 1991, JGR 96, C3, 4531-4547.

! NONE.

! ----------------------------------------------------------------------

USE cudafor
USE PARKIND_WAVE, ONLY: JWIM, JWRU, JWRB

IMPLICIT NONE

REAL(KIND=JWRB) :: AKI_ICE_CUF_HOIST_NEW
REAL(KIND=JWRB), INTENT(IN) :: G, XK, DEPTH, RHOW, CITH

! ICE PROPERTIES (assumed fixed for now)
REAL(KIND=JWRB), PARAMETER :: YMICE = 5.5E+9_JWRB ! typical value of Young modulus of sea ice
REAL(KIND=JWRB), PARAMETER :: RMUICE = 0.3_JWRB ! Poisson's ratio of sea ice
REAL(KIND=JWRB), PARAMETER :: RHOI = 922.5_JWRB ! typical value of the sea ice density

! RELATIVE ERROR LIMIT OF NEWTONS METHOD.
REAL(KIND=JWRB), PARAMETER :: EBS = 0.000001_JWRB
! MAXIMUM WAVE NUMBER
REAL(KIND=JWRB), PARAMETER :: AKI_MAX = 20.0_JWRB

REAL(KIND=JWRB) :: FICSTF, RDH
REAL(KIND=JWRB) :: OM2, AKI, AKIOLD, F, FPRIME, AKID
!$acc routine seq


IF (CITH <= 0.0_JWRB) THEN
AKI = XK
ELSE
! BENDING STIFFNESS / WATER DENSITY
FICSTF = (YMICE*CITH**3 / (12*(1 - RMUICE**2))) / RHOW

! DENSITY RATIO * ICE THICKNESS
RDH = (RHOI / RHOW)*CITH

! SQUARE OF THE OPEN OCEAN ANGULAR FREQUENCY
OM2 = G*XK*TANH(XK*DEPTH)

!* 2. ITERATION LOOP.
! ---------------

AKIOLD = 0.0_JWRB
AKI = MIN(XK, (OM2 / MAX(FICSTF, 1.0_JWRB))**0.2_JWRB)

DO WHILE (ABS(AKI - AKIOLD) > EBS*AKIOLD .and. AKI < AKI_MAX)
AKIOLD = AKI
AKID = MIN(DEPTH*AKI, 50.0_JWRB)
F = FICSTF*AKI**5 + G*AKI - OM2*(RDH*AKI + 1. / TANH(AKID))
FPRIME = 5._JWRB*FICSTF*AKI**4 + G - OM2*(RDH - DEPTH / (SINH(AKID)**2))
AKI = AKI - F / FPRIME
! in case of overshoot because it is trying to find a very large wave number
IF (AKI <= 0.0_JWRB) AKI = AKI_MAX
END DO

END IF

AKI_ICE_CUF_HOIST_NEW = AKI

END FUNCTION AKI_ICE_CUF_HOIST_NEW
END MODULE AKI_ICE_CUF_HOIST_NEW_MOD
64 changes: 64 additions & 0 deletions src/phys-scc-cuf-hoist/chnkmin.cuf_hoist_new.F90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
! (C) Copyright 1989- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!
MODULE CHNKMIN_CUF_HOIST_NEW_MOD
CONTAINS
ATTRIBUTES(DEVICE) FUNCTION CHNKMIN_CUF_HOIST_NEW (U10, ALPHA, ALPHAMIN, CHNKMIN_U)

! ----------------------------------------------------------------------

!**** *CHNKMIN* - FUNCTION TO COMPUTE THE MINMUM CHARNOCK

!* PURPOSE.
! -------


!** INTERFACE.
! ----------

! *FUNCTION* *CHNKMIN (U10)*

! METHOD.
! -------

! CHNKMIN = ALPHAMIN + (ALPHA-ALPHAMIN)*0.5_JWRB*(1.0_JWRB-TANH(U10-A))

! EXTERNALS.
! ----------

! NONE.

! REFERENCE.
! ----------

! NONE.

! ----------------------------------------------------------------------

USE cudafor
USE PARKIND_WAVE, ONLY: JWIM, JWRU, JWRB

! ----------------------------------------------------------------------

IMPLICIT NONE

REAL(KIND=JWRB) :: CHNKMIN_CUF_HOIST_NEW
REAL(KIND=JWRB), INTENT(IN) :: U10
REAL(KIND=JWRB), INTENT(IN) :: ALPHA
REAL(KIND=JWRB), INTENT(IN) :: ALPHAMIN
REAL(KIND=JWRB), INTENT(IN) :: CHNKMIN_U
!$acc routine seq

! ----------------------------------------------------------------------


CHNKMIN_CUF_HOIST_NEW = ALPHAMIN + (ALPHA - ALPHAMIN)*0.5_JWRB*(1.0_JWRB - TANH(U10 - CHNKMIN_U))


END FUNCTION CHNKMIN_CUF_HOIST_NEW
END MODULE CHNKMIN_CUF_HOIST_NEW_MOD
Loading

0 comments on commit 0f522dd

Please sign in to comment.