text:
! *********************************COPYRIGHT************************************
! (C) Crown copyright Met Office. All rights reserved.
! For further details please refer to the file LICENCE.txt
! which you should have received as part of this distribution.
! *********************************COPYRIGHT************************************
!
! This file is part of the UM Shared Library project.
!
! The UM Shared Library is free software: you can redistribute it
! and/or modify it under the terms of the Modified BSD License, as
! published by the Open Source Initiative.
!
! The UM Shared Library is distributed in the hope that it will be
! useful, but WITHOUT ANY WARRANTY; without even the implied warranty
! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! Modified BSD License for more details.
!
! You should have received a copy of the Modified BSD License
! along with the UM Shared Library.
! If not, see <http://opensource.org/licenses/BSD-3-Clause>.
!*******************************************************************************
MODULE fruit_test_shum_latlon_eq_grids_mod

USE fruit
USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT64_T, C_INT32_T, C_FLOAT, C_DOUBLE

IMPLICIT NONE

PRIVATE

PUBLIC :: fruit_test_shum_latlon_eq_grids

!------------------------------------------------------------------------------!
! We're going to use the types from the ISO_C_BINDING module, since although   !
! the REALs aren't 100% guaranteed to correspond to the sizes we want to       !
! enforce, they should be good enough on the majority of systems.              !
!                                                                              !
! Additional protection for the case that FLOAT/DOUBLE do not conform to the   !
! sizes we expect is provided via the "precision_bomb" macro-file              !
!------------------------------------------------------------------------------!
  INTEGER, PARAMETER :: int64  = C_INT64_T
  INTEGER, PARAMETER :: int32  = C_INT32_T
  INTEGER, PARAMETER :: real64 = C_DOUBLE
  INTEGER, PARAMETER :: real32 = C_FLOAT
!------------------------------------------------------------------------------!

! The number of points in the individual test arrays
  INTEGER, PARAMETER :: n = 10
! The number of distinct test cases (poles) used
  INTEGER, PARAMETER :: cases = 6

! Test tolerances
  REAL(KIND=real64), PARAMETER :: test_tol_64 = 1.0E-5_real64
  REAL(KIND=real32), PARAMETER :: test_tol_32 = 1.0E-4_real32

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

INTERFACE sample_starting_data
MODULE PROCEDURE sample_starting_data_32, sample_starting_data_64
END INTERFACE

INTERFACE sample_rotated_data
MODULE PROCEDURE sample_rotated_data_32, sample_rotated_data_64
END INTERFACE

INTERFACE sample_wind_data
MODULE PROCEDURE sample_wind_data_32, sample_wind_data_64
END INTERFACE

INTERFACE sample_rotated_wind_data
MODULE PROCEDURE sample_rotated_wind_data_32, sample_rotated_wind_data_64
END INTERFACE

CONTAINS

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

SUBROUTINE fruit_test_shum_latlon_eq_grids

USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: OUTPUT_UNIT
USE f_shum_latlon_eq_grids_version_mod, ONLY:                                  &
                                          get_shum_latlon_eq_grids_version

IMPLICIT NONE

INTEGER(KIND=int64) :: version

version = get_shum_latlon_eq_grids_version()

WRITE(OUTPUT_UNIT, "()")
WRITE(OUTPUT_UNIT, "(A,I0)")                                                   &
    "Testing shum_latlon_eq_grids at Shumlib version: ", version

CALL run_test_case(                                                            &
    test_ll_to_eq_arg64, "ll_to_eq_arg64")
CALL run_test_case(                                                            &
    test_ll_to_eq_arg32, "ll_to_eq_arg32")
CALL run_test_case(                                                            &
    test_eq_to_ll_arg64, "eq_to_ll_arg64")
CALL run_test_case(                                                            &
    test_eq_to_ll_arg32, "eq_to_ll_arg32")
CALL run_test_case(                                                            &
    test_eq_to_ll_w_arg64, "eq_to_ll_w_arg64")
CALL run_test_case(                                                            &
    test_eq_to_ll_w_arg32, "eq_to_ll_w_arg32")
CALL run_test_case(                                                            &
    test_ll_to_eq_w_arg64, "ll_to_eq_w_arg64")
CALL run_test_case(                                                            &
    test_ll_to_eq_w_arg32, "ll_to_eq_w_arg32")

END SUBROUTINE fruit_test_shum_latlon_eq_grids

! Functions used to provide sample arrays of lat/long data - the goal here
! isn't for a numerical workout but easy to identify and work with arrays.
! The range from -90 to 90 lat and 1 to 360 long is spread across 10 values.
!------------------------------------------------------------------------------!
SUBROUTINE sample_starting_data_64(latitude, longitude)
IMPLICIT NONE
REAL(KIND=real64), INTENT(OUT) :: latitude(n)
REAL(KIND=real64), INTENT(OUT) :: longitude(n)

latitude(:)  = [ -90.0_real64, -80.0_real64, -60.0_real64, -45.0_real64,       &
                 -30.0_real64,   0.0_real64,  17.0_real64,  38.0_real64,       &
                  69.0_real64,  90.0_real64  ]
longitude(:) = [   0.0_real64,   5.0_real64,  20.0_real64,  60.0_real64,       &
                  75.0_real64,  90.0_real64,  99.0_real64, 168.0_real64,       &
                 279.0_real64,   0.0_real64  ]

END SUBROUTINE sample_starting_data_64

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

SUBROUTINE sample_starting_data_32(latitude, longitude)
IMPLICIT NONE
REAL(KIND=real32), INTENT(OUT) :: latitude(n)
REAL(KIND=real32), INTENT(OUT) :: longitude(n)

REAL(KIND=real64)   :: latitude_64(n)
REAL(KIND=real64)   :: longitude_64(n)
INTEGER(KIND=int64) :: i

CALL sample_starting_data_64(latitude_64, longitude_64)

DO i=1,n
  latitude(i) = REAL(latitude_64(i), KIND=real32)
  longitude(i) = REAL(longitude_64(i), KIND=real32)
END DO

END SUBROUTINE sample_starting_data_32

! To provide decent coverage the sample rotated data is provided for 6 cases
! based on different pole positions
!------------------------------------------------------------------------------!

SUBROUTINE sample_rotated_data_64                                              &
                     (case_index, pole_lat, pole_lon, latitude_eq, longitude_eq)
IMPLICIT NONE
INTEGER(KIND=int64), INTENT(IN ) :: case_index
REAL(KIND=real64),   INTENT(OUT) :: pole_lat
REAL(KIND=real64),   INTENT(OUT) :: pole_lon
REAL(KIND=real64),   INTENT(OUT) :: latitude_eq(n)
REAL(KIND=real64),   INTENT(OUT) :: longitude_eq(n)

SELECT CASE (case_index)
  CASE (1_int64)
    ! Note; in this case the "rotated" pole is in the same position as the
    ! regular lat-lon pole, so this case should return the starting values.
    ! *Except* for the points with latitudes of -90.0 or 90.0, the longitude
    ! is effectively undefined at these points and will be returned as 0.0
    pole_lat = 90.0_real64
    pole_lon = 180.0_real64
    latitude_eq(:)  = [ -90.0_real64,                                          &
                        -80.0_real64,                                          &
                        -60.0_real64,                                          &
                        -45.0_real64,                                          &
                        -30.0_real64,                                          &
                          0.0_real64,                                          &
                         17.0_real64,                                          &
                         38.0_real64,                                          &
                         69.0_real64,                                          &
                         90.0_real64  ]

    longitude_eq(:) = [   0.0_real64,                                          &
                          5.0_real64,                                          &
                         20.0_real64,                                          &
                         60.0_real64,                                          &
                         75.0_real64,                                          &
                         90.0_real64,                                          &
                         99.0_real64,                                          &
                        168.0_real64,                                          &
                        279.0_real64,                                          &
                          0.0_real64  ]

  CASE (2_int64)
    pole_lat = 54.0_real64
    pole_lon = 130.0_real64
    latitude_eq(:)  = [ -54.0_real64,                                          &
                        -58.7895666440831_real64,                              &
                        -53.2397218048810_real64,                              &
                        -25.4617609512455_real64,                              &
                         -6.46157700979021_real64,                             &
                         26.7609842064548_real64,                              &
                         45.9183375303789_real64,                              &
                         59.6633088709791_real64,                              &
                         35.0804320498546_real64,                              &
                         54.0_real64  ]

    longitude_eq(:) = [ 179.999998090904_real64,                               &
                        164.067424780139_real64,                               &
                        128.271798984558_real64,                               &
                        132.613004167377_real64,                               &
                        134.443460342966_real64,                               &
                        133.954311410277_real64,                               &
                        134.928761741302_real64,                               &
                        253.849811448046_real64,                               &
                        346.965123607299_real64,                               &
                          0.0_real64  ]

  CASE (3_int64)
    pole_lat = 18.0_real64
    pole_lon = 80.0_real64
    latitude_eq(:)  = [ -18.0000000000000_real64,                              &
                        -15.1637492785758_real64,                              &
                         -1.71067286726960_real64,                             &
                         24.4207289101076_real64,                              &
                         41.7588127942475_real64,                              &
                         69.4894732039266_real64,                              &
                         71.8596760766394_real64,                              &
                         12.4979706108037_real64,                              &
                         -1.93508314394315_real64,                             &
                         18.0000000000000_real64  ]

    longitude_eq(:) = [ 179.999999146226_real64,                               &
                        169.992216056674_real64,                               &
                        154.328820636016_real64,                               &
                        164.596942113430_real64,                               &
                        174.192615969597_real64,                               &
                        209.709344178194_real64,                               &
                        269.751129208168_real64,                               &
                        306.230553786676_real64,                               &
                          6.70398224673472_real64,                             &
                          0.0_real64  ]

  CASE (4_int64)
    pole_lat = -18.0_real64
    pole_lon = 30.0_real64
    latitude_eq(:)  = [  18.0000000000000_real64,                              &
                         27.0005058293276_real64,                              &
                         47.3850530206082_real64,                              &
                         53.2169749313423_real64,                              &
                         47.4687959970568_real64,                              &
                         28.3937504351384_real64,                              &
                         13.6262727974178_real64,                              &
                        -48.3479264657106_real64,                              &
                        -24.2446748966316_real64,                              &
                        -18.0000000000000_real64  ]

    longitude_eq(:) = [ 179.999999146226_real64,                               &
                        175.275515243798_real64,                               &
                        172.632364760056_real64,                               &
                        216.188947557755_real64,                               &
                        244.943093060330_real64,                               &
                        280.115766071256_real64,                               &
                        293.271297230166_real64,                               &
                        307.497867496828_real64,                               &
                         21.5263019525786_real64,                              &
                          0.0_real64  ]

  CASE (5_int64)
    pole_lat = -54.0_real64
    pole_lon = -20.0_real64
    latitude_eq(:)  = [  54.0000000000000_real64,                              &
                         62.7767838116441_real64,                              &
                         67.7838555403794_real64,                              &
                         40.1082874036865_real64,                              &
                         21.1089785738208_real64,                              &
                        -11.5974542206932_real64,                              &
                        -30.6003087235067_real64,                              &
                        -73.0882699868836_real64,                              &
                        -40.7803149924612_real64,                              &
                        -54.0000000000000_real64  ]

    longitude_eq(:) = [  180.000001909096_real64,                              &
                         189.231457280872_real64,                              &
                         238.213796851124_real64,                              &
                         245.572391211510_real64,                              &
                         247.636063981486_real64,                              &
                         253.592507609032_real64,                              &
                         256.342925265174_real64,                              &
                         157.851738816216_real64,                              &
                          24.4519350174060_real64,                             &
                           0.0_real64  ]

  CASE (6_int64)
    ! Note; for the points with latitudes of 90.0 or -90.0; the longitude
    ! is effectively undefined at these points and will be returned as 0.0
    pole_lat = -90.0_real64
    pole_lon = -70.0_real64
    latitude_eq(:)  = [  90.0_real64,                                          &
                         80.0_real64,                                          &
                         60.0_real64,                                          &
                         45.0_real64,                                          &
                         30.0_real64,                                          &
                          0.0_real64,                                          &
                        -17.0_real64,                                          &
                        -38.0_real64,                                          &
                        -69.0_real64,                                          &
                        -90.0_real64  ]

    longitude_eq(:) = [  250.0_real64,                                         &
                         285.0_real64,                                         &
                         270.0_real64,                                         &
                         230.0_real64,                                         &
                         215.0_real64,                                         &
                         200.0_real64,                                         &
                         191.0_real64,                                         &
                         122.0_real64,                                         &
                          11.0_real64,                                         &
                         250.0_real64  ]

END SELECT

END SUBROUTINE sample_rotated_data_64

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

SUBROUTINE sample_rotated_data_32                                              &
                     (case_index, pole_lat, pole_lon, latitude_eq, longitude_eq)
IMPLICIT NONE
INTEGER(KIND=int32), INTENT(IN ) :: case_index
REAL(KIND=real32),   INTENT(OUT) :: pole_lat
REAL(KIND=real32),   INTENT(OUT) :: pole_lon
REAL(KIND=real32),   INTENT(OUT) :: latitude_eq(n)
REAL(KIND=real32),   INTENT(OUT) :: longitude_eq(n)

INTEGER(KIND=int64) :: case_index_64
REAL(KIND=real64)   :: pole_lat_64
REAL(KIND=real64)   :: pole_lon_64
REAL(KIND=real64)   :: latitude_eq_64(n)
REAL(KIND=real64)   :: longitude_eq_64(n)

INTEGER(KIND=int64) :: i

case_index_64 = INT(case_index, KIND=int64)

CALL sample_rotated_data_64(                                                   &
       case_index_64, pole_lat_64, pole_lon_64, latitude_eq_64, longitude_eq_64)

pole_lat = REAL(pole_lat_64, KIND=real32)
pole_lon = REAL(pole_lon_64, KIND=real32)

DO i = 1,n
  latitude_eq(i) = REAL(latitude_eq_64(i), KIND=real32)
  longitude_eq(i) = REAL(longitude_eq_64(i), KIND=real32)
END DO

END SUBROUTINE sample_rotated_data_32

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

SUBROUTINE sample_rotated_wind_data_64                                         &
                     (case_index, pole_lat, pole_lon, longitude_eq, u_ll, v_ll)
IMPLICIT NONE

INTEGER(KIND=int64), INTENT(IN ) :: case_index
REAL(KIND=real64),   INTENT(OUT) :: pole_lat
REAL(KIND=real64),   INTENT(OUT) :: pole_lon
REAL(KIND=real64),   INTENT(OUT) :: longitude_eq(n)
REAL(KIND=real64),   INTENT(OUT) :: u_ll(n)
REAL(KIND=real64),   INTENT(OUT) :: v_ll(n)

REAL(KIND=real64) :: latitude_eq_dummy(n)

! We're sharing data with the test cases from the non-vector rotation tests,
! so to avoid duplication retrieve the pole co-ordinates and longitudes
! directly from there (the latitudes aren't needed however)
CALL sample_rotated_data(                                                      &
                case_index, pole_lat, pole_lon, latitude_eq_dummy, longitude_eq)

! Now provide the expected results for the wind vectors on the eq grid
SELECT CASE (case_index)
  CASE (1_int64)
    u_ll(:) = [ -20.0_real64,                                                  &
                -12.0_real64,                                                  &
                 -5.0_real64,                                                  &
                  0.0_real64,                                                  &
                  0.3_real64,                                                  &
                  5.0_real64,                                                  &
                  9.9_real64,                                                  &
                 21.0_real64,                                                  &
                 50.0_real64,                                                  &
                120.0_real64  ]

    v_ll(:) = [ -20.0_real64,                                                  &
                 12.0_real64,                                                  &
                 -5.0_real64,                                                  &
                  0.0_real64,                                                  &
                  0.4_real64,                                                  &
                 12.0_real64,                                                  &
                 -9.9_real64,                                                  &
                 21.0_real64,                                                  &
                -50.0_real64,                                                  &
                160.0_real64  ]

  CASE (2_int64)
    u_ll(:) = [  -2.46513572980279_real64,                                     &
                 16.4036524439934_real64,                                      &
                  2.82849141954085_real64,                                     &
                  0.0_real64,                                                  &
                  0.283582564804992_real64,                                    &
                 11.7023568817882_real64,                                      &
                 -8.99897206544481_real64,                                     &
                  0.933171791845741_real64,                                    &
                 28.9585133557067_real64,                                      &
                199.701624425298_real64  ]

    v_ll(:) = [  28.1766411382487_real64,                                      &
                 -4.34973407194856_real64,                                     &
                  6.48071263747929_real64,                                     &
                  0.0_real64,                                                  &
                 -0.411802065147264_real64,                                    &
                 -5.66169969278346_real64,                                     &
                -10.7255995943883_real64,                                      &
                -29.6838203472347_real64,                                      &
                 64.5089490243591_real64,                                      &
                 10.9206777215090_real64  ]

  CASE (3_int64)
    u_ll(:) = [ -23.1691183718383_real64,                                      &
                 12.3204184796072_real64,                                      &
                 -2.89527772301524_real64,                                     &
                  0.0_real64,                                                  &
                  0.300513599726356_real64,                                    &
                 10.4673622694327_real64,                                      &
                -10.8624367341968_real64,                                      &
                  3.06913674758749_real64,                                     &
                 30.2050787224513_real64,                                      &
                136.731456986939_real64  ]

    v_ll(:) = [  16.2231918521525_real64,                                      &
                 11.6707878263360_real64,                                      &
                  6.45115237043830_real64,                                     &
                  0.0_real64,                                                  &
                 -0.399614299419373_real64,                                    &
                 -7.70936619447131_real64,                                     &
                 -8.83331495472522_real64,                                     &
                -29.5394718914304_real64,                                      &
                 63.9347575217934_real64,                                      &
               -145.960640825631_real64  ]

  CASE (4_int64)
    u_ll(:) = [ -27.3205081847729_real64,                                      &
                 -3.39028557311053_real64,                                     &
                 -6.98028101148540_real64,                                     &
                  0.0_real64,                                                  &
                  0.412019531672704_real64,                                    &
                 10.0291878814405_real64,                                      &
                -12.5568577930325_real64,                                      &
                 22.2997466637935_real64,                                      &
                 44.2103465842125_real64,                                      &
                -23.9230514129678_real64  ]

    v_ll(:) = [  -7.32050766858131_real64,                                     &
                 16.6284684722544_real64,                                      &
                 -1.12945872022677_real64,                                     &
                  0.0_real64,                                                  &
                 -0.283266513095239_real64,                                    &
                 -8.27135964873773_real64,                                     &
                 -6.19235877992844_real64,                                     &
                -19.6143136186467_real64,                                      &
                 55.1855529545896_real64,                                      &
               -198.564064249029_real64  ]

  CASE (5_int64)
    u_ll(:) = [ -25.6342556805206_real64,                                      &
                 -3.41399144132580_real64,                                     &
                 -4.58978101346331_real64,                                     &
                  0.0_real64,                                                  &
                  2.413148739746743E-002_real64,                               &
                  3.04247017090075_real64,                                     &
                -13.9246276848847_real64,                                      &
                -15.3868164870599_real64,                                      &
                  3.49840347379222_real64,                                     &
                -58.0398979394267_real64  ]

    v_ll(:) = [ -11.9534486950711_real64,                                      &
                 16.6236176098512_real64,                                      &
                  5.37902502768408_real64,                                     &
                  0.0_real64,                                                  &
                 -0.499417343748608_real64,                                    &
                -12.6389625863510_real64,                                      &
                  1.45764492629112_real64,                                     &
                -25.4016904633047_real64,                                      &
                 70.6240835206699_real64,                                      &
               -191.393234590936_real64  ]

  CASE (6_int64)
    u_ll(:) = [   2.46513666864877_real64,                                     &
                 12.0_real64,                                                  &
                  5.0_real64,                                                  &
                  0.0_real64,                                                  &
                 -0.3_real64,                                                  &
                 -5.0_real64,                                                  &
                 -9.9_real64,                                                  &
                -21.0_real64,                                                  &
                -50.0_real64,                                                  &
                 10.9206843755689_real64  ]

    v_ll(:) = [  28.1766410561103_real64,                                      &
                -12.0_real64,                                                  &
                  5.0_real64,                                                  &
                  0.0_real64,                                                  &
                 -0.4_real64,                                                  &
                -12.0_real64,                                                  &
                  9.9_real64,                                                  &
                -21.0_real64,                                                  &
                 50.0_real64,                                                  &
               -199.701624061421_real64  ]

END SELECT

END SUBROUTINE sample_rotated_wind_data_64

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

SUBROUTINE sample_rotated_wind_data_32                                         &
                     (case_index, pole_lat, pole_lon, longitude_eq, u_ll, v_ll)
IMPLICIT NONE
INTEGER(KIND=int32), INTENT(IN ) :: case_index
REAL(KIND=real32),   INTENT(OUT) :: pole_lat
REAL(KIND=real32),   INTENT(OUT) :: pole_lon
REAL(KIND=real32),   INTENT(OUT) :: longitude_eq(n)
REAL(KIND=real32),   INTENT(OUT) :: u_ll(n)
REAL(KIND=real32),   INTENT(OUT) :: v_ll(n)

INTEGER(KIND=int64) :: case_index_64
REAL(KIND=real64)   :: pole_lat_64
REAL(KIND=real64)   :: pole_lon_64
REAL(KIND=real64)   :: longitude_eq_64(n)
REAL(KIND=real64)   :: u_ll_64(n)
REAL(KIND=real64)   :: v_ll_64(n)

INTEGER(KIND=int64) :: i

case_index_64 = INT(case_index, KIND=int64)

CALL sample_rotated_wind_data_64(                                              &
     case_index_64, pole_lat_64, pole_lon_64, longitude_eq_64, u_ll_64, v_ll_64)

pole_lat = REAL(pole_lat_64, KIND=real32)
pole_lon = REAL(pole_lon_64, KIND=real32)

DO i = 1,n
  longitude_eq(i) = REAL(longitude_eq_64(i), KIND=real32)
  u_ll(i) = REAL(u_ll_64(i), KIND=real32)
  v_ll(i) = REAL(v_ll_64(i), KIND=real32)
END DO

END SUBROUTINE sample_rotated_wind_data_32

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

SUBROUTINE sample_wind_data_64(lambda, u_eq, v_eq)
IMPLICIT NONE
REAL(KIND=real64), INTENT(OUT) :: lambda(n)
REAL(KIND=real64), INTENT(OUT) :: u_eq(n)
REAL(KIND=real64), INTENT(OUT) :: v_eq(n)

REAL(KIND=real64) :: phi_dummy(n)

! We're sharing data with the test cases from the non-vector rotation tests,
! so to avoid duplication retrieve the starting longitudes directly from
! there (the latitudes aren't needed however)
CALL sample_starting_data(phi_dummy, lambda)

u_eq(:) = [ -20.0_real64, -12.0_real64,  -5.0_real64,   0.0_real64,            &
              0.3_real64,   5.0_real64,   9.9_real64,  21.0_real64,            &
             50.0_real64, 120.0_real64  ]

v_eq(:) = [ -20.0_real64,  12.0_real64,  -5.0_real64,   0.0_real64,            &
              0.4_real64,  12.0_real64,  -9.9_real64,  21.0_real64,            &
            -50.0_real64, 160.0_real64  ]

END SUBROUTINE sample_wind_data_64

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

SUBROUTINE sample_wind_data_32(lambda, u_eq, v_eq)
IMPLICIT NONE
REAL(KIND=real32), INTENT(OUT) :: lambda(n)
REAL(KIND=real32), INTENT(OUT) :: u_eq(n)
REAL(KIND=real32), INTENT(OUT) :: v_eq(n)

REAL(KIND=real64)   :: lambda64(n)
REAL(KIND=real64)   :: u_eq_64(n)
REAL(KIND=real64)   :: v_eq_64(n)

INTEGER(KIND=int64) :: i

CALL sample_wind_data_64(lambda64, u_eq_64, v_eq_64)

DO i = 1,n
  lambda(i) = REAL(lambda64(i), KIND=real32)
  u_eq(i) = REAL(u_eq_64(i), KIND=real32)
  v_eq(i) = REAL(v_eq_64(i), KIND=real32)
END DO

END SUBROUTINE sample_wind_data_32

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

SUBROUTINE test_ll_to_eq_arg64

USE f_shum_latlon_eq_grids_mod, ONLY: f_shum_latlon_to_eq

IMPLICIT NONE

REAL(KIND=real64) :: longitude(n)
REAL(KIND=real64) :: latitude(n)

REAL(KIND=real64) :: phi_pole
REAL(KIND=real64) :: phi_eq(n)
REAL(KIND=real64) :: phi_exp(n)

REAL(KIND=real64) :: lambda_pole
REAL(KIND=real64) :: lambda_eq(n)
REAL(KIND=real64) :: lambda_exp(n)

INTEGER(KIND=int64)          :: i
INTEGER(KIND=int64)          :: status
CHARACTER(LEN=500)           :: message
CHARACTER(LEN=20)            :: case_info

! Retrieve the set of lat/lon points to be tested
CALL sample_starting_data(latitude, longitude)

! The above points will be tested using a variety of target grids, the
! pole co-ordinates of these and the expected values after rotation are
! returned by a call to "sample_rotated_data" with the case number, i
DO i = 1,cases

  CALL sample_rotated_data(i, phi_pole, lambda_pole, phi_exp, lambda_exp)

  status = f_shum_latlon_to_eq(                                                &
         latitude, longitude, phi_eq, lambda_eq, phi_pole, lambda_pole, message)

  WRITE(case_info, "(A,I0,A)") " (Test pole case: ", i,")"

  CALL assert_equals(0_int64, status,                                          &
    "Rotation of ll arrays to eq returned non-zero exit status"//case_info//   &
    " Message: "//TRIM(message))

  CALL assert_equals(phi_exp, phi_eq, n, test_tol_64,                          &
    "Rotated phi array does not agree with expected result"//case_info)

  CALL assert_equals(lambda_exp, lambda_eq, n, test_tol_64,                    &
    "Rotated lambda array does not agree with expected result"//case_info)

END DO

END SUBROUTINE test_ll_to_eq_arg64

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

SUBROUTINE test_ll_to_eq_arg32

USE f_shum_latlon_eq_grids_mod, ONLY: f_shum_latlon_to_eq

IMPLICIT NONE

REAL(KIND=real32) :: longitude(n)
REAL(KIND=real32) :: latitude(n)

REAL(KIND=real32) :: phi_pole
REAL(KIND=real32) :: phi_eq(n)
REAL(KIND=real32) :: phi_exp(n)

REAL(KIND=real32) :: lambda_pole
REAL(KIND=real32) :: lambda_eq(n)
REAL(KIND=real32) :: lambda_exp(n)

INTEGER(KIND=int32)          :: i
INTEGER(KIND=int32)          :: status
CHARACTER(LEN=500)           :: message
CHARACTER(LEN=20)            :: case_info

! Retrieve the set of lat/lon points to be tested
CALL sample_starting_data(latitude, longitude)

! The above points will be tested using a variety of target grids, the
! pole co-ordinates of these and the expected values after rotation are
! returned by a call to "sample_rotated_data" with the case number, i
DO i = 1,cases

  CALL sample_rotated_data(i, phi_pole, lambda_pole, phi_exp, lambda_exp)

  status = f_shum_latlon_to_eq(                                                &
         latitude, longitude, phi_eq, lambda_eq, phi_pole, lambda_pole, message)

  WRITE(case_info, "(A,I0,A)") " (Test pole case: ", i,")"

  CALL assert_equals(0_int32, status,                                          &
    "Rotation of ll arrays to eq returned non-zero exit status"//case_info//   &
    " Message: "//TRIM(message))

  CALL assert_equals(phi_exp, phi_eq, n, test_tol_32,                          &
    "Rotated phi array does not agree with expected result"//case_info)

  CALL assert_equals(lambda_exp, lambda_eq, n, test_tol_32,                    &
    "Rotated lambda array does not agree with expected result"//case_info)

END DO

END SUBROUTINE test_ll_to_eq_arg32

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

SUBROUTINE test_eq_to_ll_arg64

USE f_shum_latlon_eq_grids_mod, ONLY: f_shum_eq_to_latlon

IMPLICIT NONE

REAL(KIND=real64) :: phi_pole
REAL(KIND=real64) :: phi_eq(n)
REAL(KIND=real64) :: phi_exp(n)

REAL(KIND=real64) :: lambda_pole
REAL(KIND=real64) :: lambda_eq(n)
REAL(KIND=real64) :: lambda_exp(n)

REAL(KIND=real64) :: latitude(n)
REAL(KIND=real64) :: longitude(n)

INTEGER(KIND=int64)          :: i
INTEGER(KIND=int64)          :: status
CHARACTER(LEN=500)           :: message
CHARACTER(LEN=20)            :: case_info

! Retrive the *original* data (which we are going to try and rotate back to)
CALL sample_starting_data(phi_exp, lambda_exp)

! Now go through the sample cases un-rotating the data and comparing to the
! original data above
DO i = 1,cases

  CALL sample_rotated_data(i, phi_pole, lambda_pole, phi_eq, lambda_eq)

  status = f_shum_eq_to_latlon(                                                &
         phi_eq, lambda_eq, latitude, longitude, phi_pole, lambda_pole, message)

  WRITE(case_info, "(A,I0,A)") " (Test pole case: ", i,")"

  CALL assert_equals(0_int64, status,                                          &
    "Unrotating eq arrays to ll returned non-zero exit status"//case_info//    &
    " Message: "//TRIM(message))

  CALL assert_equals(phi_exp, latitude, n, test_tol_64,                        &
    "Unrotated phi array does not agree with expected result"//case_info)

  CALL assert_equals(lambda_exp, longitude, n, test_tol_64,                    &
    "Unrotated lambda array does not agree with expected result"//case_info)

END DO

END SUBROUTINE test_eq_to_ll_arg64

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

SUBROUTINE test_eq_to_ll_arg32

USE f_shum_latlon_eq_grids_mod, ONLY: f_shum_eq_to_latlon

IMPLICIT NONE

REAL(KIND=real32) :: phi_pole
REAL(KIND=real32) :: phi_eq(n)
REAL(KIND=real32) :: phi_exp(n)

REAL(KIND=real32) :: lambda_pole
REAL(KIND=real32) :: lambda_eq(n)
REAL(KIND=real32) :: lambda_exp(n)

REAL(KIND=real32) :: latitude(n)
REAL(KIND=real32) :: longitude(n)

INTEGER(KIND=int32)          :: i
INTEGER(KIND=int32)          :: status
CHARACTER(LEN=500)           :: message
CHARACTER(LEN=20)            :: case_info

! Retrive the *original* data (which we are going to try and rotate back to)
CALL sample_starting_data(phi_exp, lambda_exp)

! Now go through the sample cases un-rotating the data and comparing to the
! original data above
DO i = 1,cases

  CALL sample_rotated_data(i, phi_pole, lambda_pole, phi_eq, lambda_eq)

  status = f_shum_eq_to_latlon(                                                &
         phi_eq, lambda_eq, latitude, longitude, phi_pole, lambda_pole, message)

  WRITE(case_info, "(A,I0,A)") " (Test pole case: ", i,")"

  CALL assert_equals(0_int32, status,                                          &
    "Unrotating eq arrays to ll returned non-zero exit status"//case_info//    &
    " Message: "//TRIM(message))

  CALL assert_equals(phi_exp, latitude, n, test_tol_32,                        &
    "Unrotated phi array does not agree with expected result"//case_info)

  CALL assert_equals(lambda_exp, longitude, n, test_tol_32,                    &
    "Unrotated lambda array does not agree with expected result"//case_info)

END DO

END SUBROUTINE test_eq_to_ll_arg32

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

SUBROUTINE test_eq_to_ll_w_arg64

USE f_shum_latlon_eq_grids_mod, ONLY:                                          &
                     f_shum_latlon_eq_vector_coeff, f_shum_eq_to_latlon_vector

IMPLICIT NONE

REAL(KIND=real64) :: lambda(n)
REAL(KIND=real64) :: lambda_eq(n)
REAL(KIND=real64) :: lambda_pole

REAL(KIND=real64) :: phi_pole

REAL(KIND=real64) :: coeff1(n)
REAL(KIND=real64) :: coeff2(n)

REAL(KIND=real64) :: u(n)
REAL(KIND=real64) :: u_eq(n)
REAL(KIND=real64) :: u_ll(n)

REAL(KIND=real64) :: v(n)
REAL(KIND=real64) :: v_eq(n)
REAL(KIND=real64) :: v_ll(n)

INTEGER(KIND=int64)          :: i
INTEGER(KIND=int64)          :: status
CHARACTER(LEN=500)           :: message
CHARACTER(LEN=20)            :: case_info

CALL sample_wind_data(lambda, u_eq, v_eq)

DO i = 1,cases
  CALL sample_rotated_wind_data(i, phi_pole, lambda_pole, lambda_eq, u_ll, v_ll)

  status = f_shum_latlon_eq_vector_coeff(coeff1, coeff2, lambda_eq, lambda,    &
                                         phi_pole, lambda_pole, message)

  WRITE(case_info, "(A,I0,A)") " (Test pole case: ", i,")"

  CALL assert_equals(0_int64, status,                                          &
    "Calculating vector coeffs returned non-zero exit status"//case_info//     &
    " Message: "//TRIM(message))

  status = f_shum_eq_to_latlon_vector(coeff1, coeff2, u_eq, v_eq, u, v, message)

  CALL assert_equals(0_int64, status,                                          &
       "Rotating vectors eq -> ll returned non-zero exit status"//case_info//  &
       " Message: "//TRIM(message))

  CALL assert_equals(u_ll, u, n, test_tol_64,                                  &
        "Rotated U components do not agree with expected result"//case_info)

  CALL assert_equals(v_ll, v, n, test_tol_64,                                  &
        "Rotated V components do not agree with expected result"//case_info)

END DO

END SUBROUTINE test_eq_to_ll_w_arg64

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

SUBROUTINE test_eq_to_ll_w_arg32

USE f_shum_latlon_eq_grids_mod, ONLY:                                          &
                     f_shum_latlon_eq_vector_coeff, f_shum_eq_to_latlon_vector

IMPLICIT NONE

REAL(KIND=real32) :: lambda(n)
REAL(KIND=real32) :: lambda_eq(n)
REAL(KIND=real32) :: lambda_pole

REAL(KIND=real32) :: phi_pole

REAL(KIND=real32) :: coeff1(n)
REAL(KIND=real32) :: coeff2(n)

REAL(KIND=real32) :: u(n)
REAL(KIND=real32) :: u_eq(n)
REAL(KIND=real32) :: u_ll(n)

REAL(KIND=real32) :: v(n)
REAL(KIND=real32) :: v_eq(n)
REAL(KIND=real32) :: v_ll(n)

INTEGER(KIND=int32)          :: i
INTEGER(KIND=int32)          :: status
CHARACTER(LEN=500)           :: message
CHARACTER(LEN=20)            :: case_info

CALL sample_wind_data(lambda, u_eq, v_eq)

DO i = 1,cases
  CALL sample_rotated_wind_data(i, phi_pole, lambda_pole, lambda_eq, u_ll, v_ll)

  status = f_shum_latlon_eq_vector_coeff(coeff1, coeff2, lambda_eq, lambda,    &
                                         phi_pole, lambda_pole, message)

  WRITE(case_info, "(A,I0,A)") " (Test pole case: ", i,")"

  CALL assert_equals(0_int32, status,                                          &
    "Calculating vector coeffs returned non-zero exit status"//case_info//     &
    " Message:"//TRIM(message))

  status = f_shum_eq_to_latlon_vector(coeff1, coeff2, u_eq, v_eq, u, v, message)

  CALL assert_equals(0_int32, status,                                          &
       "Rotating vectors eq -> ll returned non-zero exit status"//case_info//  &
       " Message:"//TRIM(message))

  CALL assert_equals(u_ll, u, n, test_tol_32,                                  &
        "Rotated U components do not agree with expected result"//case_info)

  CALL assert_equals(v_ll, v, n, test_tol_32,                                  &
        "Rotated V components do not agree with expected result"//case_info)

END DO

END SUBROUTINE test_eq_to_ll_w_arg32

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

SUBROUTINE test_ll_to_eq_w_arg64

USE f_shum_latlon_eq_grids_mod, ONLY:                                          &
                     f_shum_latlon_eq_vector_coeff, f_shum_latlon_to_eq_vector

IMPLICIT NONE

REAL(KIND=real64) :: lambda(n)
REAL(KIND=real64) :: lambda_eq(n)
REAL(KIND=real64) :: lambda_pole

REAL(KIND=real64) :: phi_pole

REAL(KIND=real64) :: coeff1(n)
REAL(KIND=real64) :: coeff2(n)

REAL(KIND=real64) :: u(n)
REAL(KIND=real64) :: u_eq(n)
REAL(KIND=real64) :: u_ll(n)

REAL(KIND=real64) :: v(n)
REAL(KIND=real64) :: v_eq(n)
REAL(KIND=real64) :: v_ll(n)

INTEGER(KIND=int64)          :: i
INTEGER(KIND=int64)          :: status
CHARACTER(LEN=500)           :: message
CHARACTER(LEN=20)            :: case_info

CALL sample_wind_data(lambda, u_eq, v_eq)

DO i = 1,cases
  CALL sample_rotated_wind_data(i, phi_pole, lambda_pole, lambda_eq, u_ll, v_ll)

  status = f_shum_latlon_eq_vector_coeff(coeff1, coeff2, lambda_eq, lambda,    &
                                         phi_pole, lambda_pole, message)

  WRITE(case_info, "(A,I0,A)") " (Test pole case: ", i,")"

  CALL assert_equals(0_int64, status,                                          &
      "Calculating vector coeffs returned non-zero exit status"//case_info//   &
      " Message: "//TRIM(message))

  status = f_shum_latlon_to_eq_vector(coeff1, coeff2, u_ll, v_ll, u, v, message)

  CALL assert_equals(0_int64, status,                                          &
      "Rotating vectors ll -> eq returned non-zero exit status"//case_info//   &
      " Message: "//TRIM(message))

  CALL assert_equals(u_eq, u, n, test_tol_64,                                  &
      "Unrotated U components do not agree with expected result"//case_info)

  CALL assert_equals(v_eq, v, n, test_tol_64,                                  &
      "Unrotated V components do not agree with expected result"//case_info)

END DO

END SUBROUTINE test_ll_to_eq_w_arg64

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

SUBROUTINE test_ll_to_eq_w_arg32

USE f_shum_latlon_eq_grids_mod, ONLY:                                          &
                     f_shum_latlon_eq_vector_coeff, f_shum_latlon_to_eq_vector

IMPLICIT NONE

REAL(KIND=real32) :: lambda(n)
REAL(KIND=real32) :: lambda_eq(n)
REAL(KIND=real32) :: lambda_pole

REAL(KIND=real32) :: phi_pole

REAL(KIND=real32) :: coeff1(n)
REAL(KIND=real32) :: coeff2(n)

REAL(KIND=real32) :: u(n)
REAL(KIND=real32) :: u_eq(n)
REAL(KIND=real32) :: u_ll(n)

REAL(KIND=real32) :: v(n)
REAL(KIND=real32) :: v_eq(n)
REAL(KIND=real32) :: v_ll(n)

INTEGER(KIND=int32)          :: i
INTEGER(KIND=int32)          :: status
CHARACTER(LEN=500)           :: message
CHARACTER(LEN=20)            :: case_info

CALL sample_wind_data(lambda, u_eq, v_eq)

DO i = 1,cases
  CALL sample_rotated_wind_data(i, phi_pole, lambda_pole, lambda_eq, u_ll, v_ll)

  status = f_shum_latlon_eq_vector_coeff(coeff1, coeff2, lambda_eq, lambda,    &
                                         phi_pole, lambda_pole, message)

  WRITE(case_info, "(A,I0,A)") " (Test pole case: ", i,")"

  CALL assert_equals(0_int32, status,                                          &
      "Calculating vector coeffs returned non-zero exit status"//case_info//   &
      " Message: "//TRIM(message))

  status = f_shum_latlon_to_eq_vector(coeff1, coeff2, u_ll, v_ll, u, v, message)

  CALL assert_equals(0_int32, status,                                          &
      "Rotating vectors ll -> eq returned non-zero exit status"//case_info//   &
      " Message: "//TRIM(message))

  CALL assert_equals(u_eq, u, n, test_tol_32,                                  &
      "Unrotated U components do not agree with expected result"//case_info)

  CALL assert_equals(v_eq, v, n, test_tol_32,                                  &
      "Unrotated V components do not agree with expected result"//case_info)

END DO

END SUBROUTINE test_ll_to_eq_w_arg32

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

END MODULE fruit_test_shum_latlon_eq_grids_mod

max_stars_repo_name: climbfuji/shumlib

max_issues_repo_name: climbfuji/shumlib

max_forks_repo_name: climbfuji/shumlib


text:
! Creates a new local Git repository.
program main
    use, intrinsic :: iso_c_binding
    use :: git2
    implicit none
    character(len=*), parameter :: GIT_DIR = '.'

    integer     :: rc
    type(c_ptr) :: repo

    print '(a)', 'Init ...'
    rc = git_libgit2_init()
    if (rc < 0) call print_last_error()

    print '(3a)', 'Init repository "', GIT_DIR, '" ...'
    rc = git_repository_init(repo, GIT_DIR // c_null_char, .false._c_bool)
    if (rc < 0) call print_last_error()

    print '(a)', 'Shutdown ...'
    rc = git_libgit2_shutdown()
    if (rc < 0) call print_last_error()
contains
    subroutine print_last_error()
        character(len=:), allocatable :: message
        type(git_error), pointer      :: error

        error => git_error_last()
        call c_f_str_ptr(error%message, message)
        write (*, '("Error: ", a)') message
    end subroutine print_last_error
end program main

max_stars_repo_name: interkosmos/fortran-git

max_issues_repo_name: interkosmos/fortran-git

max_forks_repo_name: interkosmos/fortran-git


text:
!==============================================================================!
  subroutine Load_Domain(dom, grid)
!------------------------------------------------------------------------------!
!   Reads: .dom file                                                           !
!----------------------------------[Modules]-----------------------------------!
  use all_mod
  use gen_mod
  use Tokenizer_Mod
  use Domain_Mod
  use Grid_Mod
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Domain_Type) :: dom
  type(Grid_Type)   :: grid
!----------------------------------[Calling]-----------------------------------!
  real :: Tet_Volume
!-----------------------------------[Locals]-----------------------------------!
  integer           :: b, i, l, s, fc, n, n1,n2,n3,n4
  integer           :: n_faces_check, n_nodes_check
  integer           :: ni, nj, nk, npnt
  character(len=12) :: dum
  character(len=80) :: domain_name
  character(len=12) :: answer
  real              :: xt(8), yt(8), zt(8)
  integer           :: face_nodes(6,4)
!==============================================================================!
  data face_nodes / 1, 1, 2, 4, 3, 5,         &
                    2, 5, 6, 8, 7, 7,         &
                    4, 6, 8, 7, 5, 8,         &
                    3, 2, 4, 3, 1, 6  /
!------------------------------------------------------------------------------!

  print *, '#========================================'
  print *, '# Input problem name: (without extension)'
  print *, '#----------------------------------------'
  call Tokenizer_Mod_Read_Line(5)
  read(line % tokens(1), *) name

  domain_name = name
  domain_name(len_trim(name)+1:len_trim(name)+4) = '.dom'
  write(*, *) '# Now reading the file: ', domain_name
  open(9, file=domain_name)

  !-----------------------------------------------------------------!
  !   Max. number of nodes (cells), boundary faces and cell faces   !
  !-----------------------------------------------------------------!
  call Tokenizer_Mod_Read_Line(9)
  read(line % tokens(1), *) grid % max_n_nodes
  read(line % tokens(2), *) grid % max_n_bnd_cells
  read(line % tokens(3), *) grid % max_n_faces

  !---------------------!
  !   Allocate memory   !
  !---------------------!
  print *, '# Allocating memory for: '
  print *, '#', grid % max_n_nodes, ' nodes and cells'
  print *, '#', grid % max_n_bnd_cells, ' boundary cells'
  print *, '#', grid % max_n_faces, ' cell faces'

  ! Variables declared in all_mod.h90:
  allocate (delta (-grid % max_n_bnd_cells:grid % max_n_nodes))
  delta=0.0

  allocate (WallDs(grid % max_n_nodes))
  WallDs=0.0

  allocate (material(-grid % max_n_bnd_cells:grid % max_n_nodes))
  material=0

  allocate (CopyC(-grid % max_n_bnd_cells:grid % max_n_nodes))
  CopyC=0
  allocate (CopyS(2,grid % max_n_bnd_cells))
  CopyS=0

  allocate (BCmark(-grid % max_n_bnd_cells-1:-1))
  BCmark=0;

  ! Variables in Grid_Mod
  call Grid_Mod_Allocate_Nodes(grid,  &
                               grid % max_n_nodes)

  call Grid_Mod_Allocate_Cells(grid,                    &
                               grid % max_n_nodes,      &
                               grid % max_n_bnd_cells)

  call Grid_Mod_Allocate_Faces(grid,                &
                               grid % max_n_faces,  &
                               0)                      ! no shadow faces

  ! Variables declared in gen_mod.h90:
  allocate (walln(grid % max_n_nodes))
  walln=0

  allocate (SideCc(grid % max_n_faces,2))
  SideCc=0

  ! Variables for renumbering
  allocate (NewN(-grid % max_n_bnd_cells:grid % max_n_nodes))
  NewN=0
  allocate (NewC(-grid % max_n_bnd_cells:grid % max_n_nodes))
  NewC=0
  allocate (NewS( grid % max_n_faces))
  NewS=0
  allocate (CelMar(-grid % max_n_bnd_cells:grid % max_n_nodes))
  CelMar=0

  allocate (TwinN(grid % max_n_nodes,0:8))
  TwinN=0

  allocate (level(grid % max_n_nodes))
  level=0

  print *, '# Allocation successfull !'

  !-------------!
  !   Corners   !
  !-------------!
  call Tokenizer_Mod_Read_Line(9)
  read(line % tokens(1), *) dom % n_points  ! number of points

  call Domain_Mod_Allocate_Points(dom, dom % n_points)

  do i = 1, dom % n_points
    call Tokenizer_Mod_Read_Line(9)
    read(line % tokens(2),*) dom % points(i) % x
    read(line % tokens(3),*) dom % points(i) % y
    read(line % tokens(4),*) dom % points(i) % z
  end do

  !------------!
  !   Blocks   !
  !------------!
  call Tokenizer_Mod_Read_Line(9)
  read(line % tokens(1), *) dom % n_blocks  ! number of blocks

  call Domain_Mod_Allocate_Blocks(dom, dom % n_blocks)

  ! Initialize weights
  do b=1, dom % n_blocks
    dom % blocks(b) % weights      = 1.0
    dom % blocks(b) % face_weights = 1.0
  end do

  do b = 1, dom % n_blocks
    dom % blocks(b) % corners(0)=1       ! suppose it is properly oriented

    call Tokenizer_Mod_Read_Line(9)
    read(line % tokens(2),*) dom % blocks(b) % resolutions(1)
    read(line % tokens(3),*) dom % blocks(b) % resolutions(2)
    read(line % tokens(4),*) dom % blocks(b) % resolutions(3)

    call Tokenizer_Mod_Read_Line(9)
    read(line % whole, *)               &  ! block weights
         dom % blocks(b) % weights(1),  &
         dom % blocks(b) % weights(2),  &
         dom % blocks(b) % weights(3)

    call Tokenizer_Mod_Read_Line(9)
    read(line % whole, *)                                             &
         dom % blocks(b) % corners(1), dom % blocks(b) % corners(2),  &
         dom % blocks(b) % corners(3), dom % blocks(b) % corners(4),  &
         dom % blocks(b) % corners(5), dom % blocks(b) % corners(6),  &
         dom % blocks(b) % corners(7), dom % blocks(b) % corners(8)

    !---------------------------!
    !   Check if the block is   !
    !     properly oriented     !
    !---------------------------!
    do n=1,8
      xt(n) = dom % points(dom % blocks(b) % corners(n)) % x
      yt(n) = dom % points(dom % blocks(b) % corners(n)) % y
      zt(n) = dom % points(dom % blocks(b) % corners(n)) % z
    end do

    if(Tet_Volume( xt(2),yt(2),zt(2), xt(5),yt(5),zt(5),  &
                   xt(3),yt(3),zt(3), xt(1),yt(1),zt(1) )  < 0) then
      dom % blocks(b) % corners(0)=-1            !  It's nor properly oriented
      call Swap_Integers(dom % blocks(b) % corners(2),  &
                         dom % blocks(b) % corners(3))
      call Swap_Integers(dom % blocks(b) % corners(6),  &
                         dom % blocks(b) % corners(7))
      call Swap_Reals(dom % blocks(b) % weights(1),  &
                      dom % blocks(b) % weights(2))
      dom % blocks(b) % weights(1) = 1.0 / dom % blocks(b) % weights(1)
      dom % blocks(b) % weights(2) = 1.0 / dom % blocks(b) % weights(2)
      call Swap_Integers(dom % blocks(b) % resolutions(1),  &
                         dom % blocks(b) % resolutions(2))
      print *, 'Warning: Block ',b,' was not properly oriented'
    end if
  end do                 ! through dom % blocks

  !-----------------------------!
  !   Set the corners of each   !
  !      face of the block      !
  !-----------------------------!
  do b = 1, dom % n_blocks
    do fc = 1,6
      do n = 1,4
        dom % blocks(b) % faces(fc, n) =  &
        dom % blocks(b) % corners(face_nodes(fc,n))
      end do
    end do
  end do

  !----------------------------------------------!
  !   Lines                                      !
  !----------------------------------------------!
  !   Lines can be prescribed point by point     !
  !   or with just a weighting factor.           !
  !----------------------------------------------!
  call Tokenizer_Mod_Read_Line(9)
  read(line % tokens(1), *) dom % n_lines  ! number of defined dom % lines

  call Domain_Mod_Allocate_Lines(dom, dom % n_lines)

  do l=1, dom % n_lines
    call Tokenizer_Mod_Read_Line(9)

    read(line % tokens(1),*) npnt
    read(line % tokens(2),*) dom % lines(l) % points(1)
    read(line % tokens(3),*) dom % lines(l) % points(2)

    call Find_Line(dom,                         &
                   dom % lines(l) % points(1),  &
                   dom % lines(l) % points(2),  &
                   dom % lines(l) % resolution)

    ! Does this need a more elegant solution?
    allocate(dom % lines(l) % x( dom % lines(l) % resolution ))
    allocate(dom % lines(l) % y( dom % lines(l) % resolution ))
    allocate(dom % lines(l) % z( dom % lines(l) % resolution ))

    ! Point by point
    if(npnt > 0) then
      do n=1,dom % lines(l) % resolution
        call Tokenizer_Mod_Read_Line(9)
        read(line % tokens(2),*) dom % lines(l) % x(n)
        read(line % tokens(3),*) dom % lines(l) % y(n)
        read(line % tokens(4),*) dom % lines(l) % z(n)
      end do

    ! Weight factor
    else
      call Tokenizer_Mod_Read_Line(9)
      read(line % tokens(1), *) dom % lines(l) % weight
    endif

  end do

  !----------------------------------------!
  !   Copy block weights to face weights   !
  !----------------------------------------!
  do b = 1, dom % n_blocks
    do fc = 1,6                          !  face of the block
      dom % blocks(b) % face_weights(fc, 1) = dom % blocks(b) % weights(1)
      dom % blocks(b) % face_weights(fc, 2) = dom % blocks(b) % weights(2)
      dom % blocks(b) % face_weights(fc, 3) = dom % blocks(b) % weights(3)
    end do
  end do

  !--------------!
  !   Surfaces   !
  !--------------!
  call Tokenizer_Mod_Read_Line(9)
  read(line % tokens(1), *) Nsurf     ! number of defined surfaces

  do s=1,Nsurf
    call Tokenizer_Mod_Read_Line(9)
    read(line % whole,*) dum, n1, n2, n3, n4
    call Find_Surface(dom, n1, n2, n3, n4, b, fc)
    print *, '# block: ', b, ' surf: ', fc
    n = (b-1)*6 + fc         ! surface number

    call Tokenizer_Mod_Read_Line(9)
    read(line % whole, *) dom % blocks(b) % face_weights(fc,1),  &
                          dom % blocks(b) % face_weights(fc,2),  &
                          dom % blocks(b) % face_weights(fc,2)
  end do

  !---------------------------------------!
  !   Is there enough allocated memory?   !
  !---------------------------------------!

  ! Nodes & faces
  n_nodes_check = 0
  n_faces_check = 0
  do b = 1, dom % n_blocks
    ni = dom % blocks(b) % resolutions(1)
    nj = dom % blocks(b) % resolutions(2)
    nk = dom % blocks(b) % resolutions(3)
    n_nodes_check=n_nodes_check + ni*nj*nk
    n_faces_check=n_faces_check + ni*nj*nk + 2*( (ni*nj)+(nj*nk)+(ni*nk) )
  end do

  if( (n_faces_check > grid % max_n_faces) .or.  &
      (n_nodes_check > grid % max_n_nodes) ) then
    print *, '# Error message from TFlowS:'
  end if

  if( n_faces_check  > grid % max_n_faces ) then
    print *, '# The estimated number of sides is :', n_faces_check
    print *, '# There is space available only for:', grid % max_n_faces
    print *, '# Increase the parameter grid % max_n_faces in the input file'
    print *, '# and re-run the code !'
  end if

  if( n_nodes_check  > grid % max_n_nodes ) then
    print *, '# The estimated number of nodes is :', n_nodes_check
    print *, '# There is space available only for:', grid % max_n_nodes
    print *, '# Increase the parameter grid % max_n_nodes in the input file'
    print *, '# and re-run the code !'
  end if

  if( (n_faces_check > grid % max_n_faces) .or.  &
      (n_nodes_check > grid % max_n_nodes) ) then
    stop
  end if

  !---------------------------------------!
  !   Boundary conditions and materials   !
  !---------------------------------------!
  call Tokenizer_Mod_Read_Line(9)
  read(line % tokens(1), *) dom % n_regions  ! number of regions (can be boundary
                                             ! conditions or materials)

  call Domain_Mod_Allocate_Regions(dom, dom % n_regions)

  do n=1, dom % n_regions
    dom % regions(n) % face=''

    call Tokenizer_Mod_Read_Line(9)
    if(line % n_tokens == 7) then
      read(line % whole,*)  dum,            &
                   dom % regions(n) % is,   &
                   dom % regions(n) % js,   &
                   dom % regions(n) % ks,   &
                   dom % regions(n) % ie,   &
                   dom % regions(n) % je,   &
                   dom % regions(n) % ke
    else if(line % n_tokens == 2) then
      read(line % tokens(1),*)       dum
      read(line % tokens(2),'(A4)')  &
           dom % regions(n) % face
      call To_Upper_Case(dom % regions(n) % face)
    end if

    call Tokenizer_Mod_Read_Line(9)
    read(line % tokens(1), *) dom % regions(n) % block
    read(line % tokens(2), *) dom % regions(n) % name
    call To_Upper_Case(dom % regions(n) % name)

    ! if( dom % blocks(b_cond(n,7)) % points(0) == -1 ) then
    !   call Swap_Integers( dom % regions(n) % is,dom % regions(n) % js )
    !   call Swap_Integers( dom % regions(n) % ie,dom % regions(n) % je )
    ! end if

  end do

  !-------------------------!
  !   Periodic boundaries   !
  !-------------------------!
  call Tokenizer_Mod_Read_Line(9)
  read(line % tokens(1), *)  n_periodic_cond      ! number of periodic boundaries
  print *, '# Number of periodic boundaries: ', n_periodic_cond

  allocate (periodic_cond(n_periodic_cond,8))

  do n=1,n_periodic_cond
    call Tokenizer_Mod_Read_Line(9)
    read(line % whole, *) dum, periodic_cond(n,1), periodic_cond(n,2),  &
                               periodic_cond(n,3), periodic_cond(n,4)
    call Tokenizer_Mod_Read_Line(9)
    read(line % whole, *)      periodic_cond(n,5), periodic_cond(n,6),  &
                               periodic_cond(n,7), periodic_cond(n,8)
  end do

  !---------------------!
  !   Copy boundaries   !
  !---------------------!
  call Tokenizer_Mod_Read_Line(9)
  read(line % tokens(1), *)  n_copy_cond  ! number of copy boundaries
  print *, '# Number of copy boundaries: ', n_copy_cond

  allocate (copy_cond(n_copy_cond,8))

  do n=1,n_copy_cond
    call Tokenizer_Mod_Read_Line(9)
    read(line % whole, *) dum, copy_cond(n,1), copy_cond(n,2),  &
                               copy_cond(n,3), copy_cond(n,4)
    call Tokenizer_Mod_Read_Line(9)
    read(line % whole, *)      copy_cond(n,5), copy_cond(n,6),  &
                               copy_cond(n,7),copy_cond(n,8)
    call Tokenizer_Mod_Read_Line(9)
    read(line % whole, *)  copy_cond(n,0)
  end do

  !-----------------------------------!
  !   Refinement levels and regions   !
  !-----------------------------------!
  call Tokenizer_Mod_Read_Line(9)
  read(line % tokens(1), *) n_refine_levels ! number of refinement levels

  print *, '# Number of refinement levels: ', n_refine_levels

  allocate (refined_regions(n_refine_levels, 1024, 0:6))
  allocate (n_refined_regions(n_refine_levels))

  do l=1,n_refine_levels
    print *, 'Level: ', l
    call Tokenizer_Mod_Read_Line(9)
    read(line % tokens(2), *) n_refined_regions(l)

    ! Browse through regions in level "l"
    do n = 1, n_refined_regions(l)
      call Tokenizer_Mod_Read_Line(9)
      read(line % tokens(3),*) answer
      call To_Upper_Case(answer)
      if(answer == 'RECTANGLE') then
        refined_regions(l,n,0) = RECTANGLE
      elseif(answer == 'ELIPSOID') then
        refined_regions(l,n,0) = ELIPSOID
      elseif(answer == 'PLANE') then
        refined_regions(l,n,0) = PLANE
      else
        print *, 'Error in input file: ', answer
        stop
      endif

      call Tokenizer_Mod_Read_Line(9)
      read(line % whole, *)                                       &
                refined_regions(l,n,1),refined_regions(l,n,2),    &
                refined_regions(l,n,3),refined_regions(l,n,4),    &
                refined_regions(l,n,5),refined_regions(l,n,6)
    end do
  end do

  !-----------------------!
  !   Smoothing regions   !
  !-----------------------!
  call Tokenizer_Mod_Read_Line(9)
  read(line % tokens(1), *) n_smoothing_regions  ! number of smoothing regions

  print *, '# Number of (non)smoothing regions: ', n_smoothing_regions

  allocate (smooth_in_x   (n_smoothing_regions))
  allocate (smooth_in_y   (n_smoothing_regions))
  allocate (smooth_in_z   (n_smoothing_regions))
  allocate (smooth_iters  (n_smoothing_regions))
  allocate (smooth_relax  (n_smoothing_regions))
  allocate (smooth_regions(n_smoothing_regions, 0:6))

  do n=1, n_smoothing_regions
    smooth_in_x(n) = .FALSE.
    smooth_in_y(n) = .FALSE.
    smooth_in_z(n) = .FALSE.
    call Tokenizer_Mod_Read_Line(9)
    read(line % tokens(1), *) smooth_regions(n,0)
    if(line % n_tokens == 4) then   ! smoothing in three directions
      smooth_in_x(n) = .TRUE.
      smooth_in_y(n) = .TRUE.
      smooth_in_z(n) = .TRUE.
    else if(line % n_tokens == 3) then
      call To_Upper_Case(line % tokens(2))
      call To_Upper_Case(line % tokens(3))
      if( line % tokens(2)  ==  'X' )  &
          smooth_in_x(n) = .TRUE.
      if( line % tokens(3)  ==  'X' )  &
          smooth_in_x(n) = .TRUE.
      if( line % tokens(2)  ==  'Y' )  &
          smooth_in_y(n) = .TRUE.
      if( line % tokens(3)  ==  'Y' )  &
          smooth_in_y(n) = .TRUE.
      if( line % tokens(2)  ==  'Z' )  &
          smooth_in_z(n) = .TRUE.
      if( line % tokens(3)  ==  'Z' )  &
          smooth_in_z(n) = .TRUE.
    else if(line % n_tokens == 2) then
      call To_Upper_Case(line % tokens(2))
      if( line % tokens(2)  ==  'X' )  &
          smooth_in_x(n) = .TRUE.
      if( line % tokens(2)  ==  'Y' )  &
          smooth_in_y(n) = .TRUE.
      if( line % tokens(2)  ==  'Z' )  &
          smooth_in_z(n) = .TRUE.
    end if

    ! Read the coordinates of the (non)smoothed region
    call Tokenizer_Mod_Read_Line(9)
    read(line % whole, *) smooth_iters(n), smooth_relax(n)
    call Tokenizer_Mod_Read_Line(9)
    read(line % whole, *) smooth_regions(n,1),  &
                          smooth_regions(n,2),  &
                          smooth_regions(n,3),  &
                          smooth_regions(n,4),  &
                          smooth_regions(n,5),  &
                          smooth_regions(n,6)
  end do

  close(9)

  end subroutine

max_stars_repo_name: mhadziabdic/T-FlowS

max_issues_repo_name: mhadziabdic/T-FlowS

max_forks_repo_name: mhadziabdic/T-FlowS


!!$ This module sets up the floating point precision used.

module real_type_mod

!!$ sp contains the kind value for single precision.
  integer, parameter :: sp = kind(1.0)

!!$ dp contains the kind value for double precision.
  integer, parameter :: dp = kind(1.0d0)

!!$ Set wp to the desired precision.
  integer, parameter :: wp = dp

end module real_type_mod


text:
subroutine calseddf2004(ustarc    ,ws        ,tp        ,hrms      ,h1        , &
                      & seddif    ,kmax      ,sig       ,thick     ,dicww     , &
                      & tauwav    ,tauc      ,ltur      ,delw      ,rhowat    , &
                      & uwbih     ,aks       ,caks      ,caks_ss3d ,deltas    , &
                      & aks_ss3d  ,d50       ,salinity  ,ws0       ,psi       , &
                      & epspar    ,eps       ,vonkar    ,salmax    ,wave      )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: calseddf2004.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/utils_gpl/morphology/packages/morphology_kernel/src/calseddf2004.f90 $
!!--description-----------------------------------------------------------------
!
! Compute sediment diffusion coefficient
! Van Rijn (2004)
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use mathconsts
    use sediment_basics_module
    !
    implicit none
!
! Call variables
!
    integer                    , intent(in)  :: kmax   !  Description and declaration in esm_alloc_int.f90
    integer                    , intent(in)  :: ltur   !  Description and declaration in esm_alloc_int.f90
    real(fp)                   , intent(in)  :: psi
    real(fp)                   , intent(in)  :: aks
    real(fp)                   , intent(out) :: aks_ss3d
    real(fp)                   , intent(in)  :: caks
    real(fp)                   , intent(out) :: caks_ss3d
    real(fp)                   , intent(in)  :: d50
    real(fp)                   , intent(in)  :: delw
    real(fp)                   , intent(out) :: deltas
    real(fp)                   , intent(in)  :: h1
    real(fp)                   , intent(in)  :: hrms   !  Description and declaration in esm_alloc_real.f90
    real(fp)                   , intent(in)  :: rhowat !  Description and declaration in esm_alloc_real.f90
    real(fp)                   , intent(in)  :: salinity
    real(fp)                   , intent(in)  :: salmax
    real(fp)                   , intent(in)  :: tauc
    real(fp)                   , intent(in)  :: tauwav
    real(fp)                   , intent(in)  :: tp     !  Description and declaration in esm_alloc_real.f90
    real(fp)                   , intent(in)  :: ustarc
    real(fp)                   , intent(in)  :: uwbih
    real(fp)                   , intent(in)  :: ws0
    real(fp), dimension(0:kmax), intent(in)  :: dicww  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(0:kmax), intent(out) :: seddif !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(0:kmax), intent(in)  :: ws     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(kmax)  , intent(in)  :: sig    !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(kmax)  , intent(in)  :: thick  !  Description and declaration in esm_alloc_real.f90
    real(fp)                   , intent(in)  :: eps
    real(fp)                   , intent(in)  :: vonkar
    logical                    , intent(in)  :: epspar
    logical                    , intent(in)  :: wave
!
! Local variables
!
    integer                     :: i
    real(fp)                    :: betaw
    real(fp)                    :: c
    real(fp)                    :: cnew
    real(fp)                    :: cmax
    real(fp)                    :: cmaxs
    real(fp)                    :: efloc
    real(fp)                    :: epsbed
    real(fp)                    :: epscur
    real(fp)                    :: epsmax
    real(fp)                    :: epsmxc
    real(fp)                    :: epstot
    real(fp)                    :: epswav
    real(fp)                    :: fch2
    real(fp)                    :: fdamp
    real(fp)                    :: ffloc
    real(fp)                    :: hinset
    real(fp)                    :: fhulp
    real(fp)                    :: gambr
    real(fp)                    :: ff
    real(fp)                    :: fi
    real(fp)                    :: hs
    real(fp)                    :: ustw
    real(fp)                    :: z
    real(fp)                    :: zkmax
    real(fp)                    :: znew
    logical                     :: difvr
!
!! executable statements -------------------------------------------------------
!
    if (tp>0.0_fp .and. wave) then
       hs = hrms * sqrt(2.0_fp)
       if (hs/h1 > 0.4_fp) then
          gambr = 1.0_fp + sqrt((hs/h1)-0.4_fp)
       else
          gambr = 1.0_fp
       endif
       deltas = 2.0_fp*gambr*delw
       deltas = min(max(0.05_fp, deltas), 0.2_fp)
       ustw   = sqrt(tauwav/rhowat)
       betaw  = min(1.0_fp+2.0_fp*(ws(1)/max(1.0e-5_fp , ustw))**2, 1.5_fp)
       !
       ! epsbed is now based on isobe-horikawa i.s.o. uwb
       ! including a damping factor as a function of the mobility parameter
       !
       fdamp  = min(max(0.1_fp , sqrt(250.0_fp/psi)), 1.0_fp)
       epsbed = 0.018_fp * fdamp * betaw * gambr * deltas * uwbih
       epsmax = 0.035_fp * gambr * h1 * hs / tp
       epsmax = min(max(epsbed , epsmax), 0.05_fp)
    else
       deltas = 0.05_fp
       epsbed = 0.0_fp
       epsmax = 0.0_fp
    endif
    !
    ! calculate vertical sediment diffusion coefficient
    !
    difvr = epspar .and. wave
    call calseddf1993(ustarc    ,ws        ,h1        ,kmax      ,sig       , &
                    & thick     ,dicww     ,tauwav    ,tauc      ,ltur      , &
                    & eps       ,vonkar    ,difvr     ,deltas    ,epsbed    , &
                    & epsmax    ,epsmxc    ,seddif    )
    !
    ! Determine height cell centre kmax
    !
    zkmax = (1.0_fp+sig(kmax)) * h1
    if (zkmax>aks .and. caks>1.0e-5_fp) then
       !
       ! If aks is lower than cell centre of kmax layer, compute new reference
       ! concentration in kmax centre by working upward from aks. This new
       ! reference concentration will be used in the source and sink terms
       ! but it will not be used in determining the suspended load correction
       ! in bott3d; for the latter the original reference height and
       ! concentration are used.
       !
       cmaxs = 0.65_fp
       cmax  = min(max(0.05_fp, (d50/dsand)*cmaxs) , cmaxs)
       ffloc = 1.0_fp
       fch2  = min(max(0.3_fp, d50/(1.5_fp*dsand)) , 1.0_fp)
       !
       c = caks
       z = aks
       do i = 1, 20
          !
          ! epstot: diffusion term
          ! note: z will always be less than or equal to h1/2
          !
          if (z > deltas) then
             epswav = epsbed + (epsmax-epsbed) * ((z-deltas)/(0.5_fp*h1-deltas))
          else
             epswav = epsbed
          endif
          epscur = epsmxc * (1.0_fp-(1.0_fp-2.0_fp*z/h1)**2)
          epstot = sqrt(epswav**2 + epscur**2)
          !
          ! fi : damping term
          !
          fi  = max(0.01_fp, fch2*(1.0_fp+((c/cmaxs)**0.8_fp)-2.0_fp*((c/cmaxs)**0.4_fp)))
          !
          ! hinset : hindered settling term
          ! note: the effect of other sediment fractions is missing, compare FALLVE
          !
          hinset = max(0.01_fp, sqrt(1.0_fp-0.65_fp*c/cmax))
          !
          ! ffloc : flocculation term
          ! note: the effect of other sediment fractions is missing, compare FALLVE
          !
          if (d50 < dsand .and. salinity>0.01_fp .and. salmax>0.0_fp) then
             fhulp = max(1.0_fp, 4.0_fp+log10(2.0_fp*c/cmax))
             efloc = min(max(0.0_fp, dsand/d50-1.0_fp) , 3.0_fp)
             ffloc = min(max(1.0_fp, fhulp**efloc) , 10.0_fp)
             ffloc = (ffloc-1.0_fp) * min(1.0_fp , salinity/salmax) + 1.0_fp
          endif
          ff   = -ws0 * hinset * ffloc / (epstot * fi)
          !
          znew = min(aks*(zkmax/aks)**(real(i,fp)/20.0_fp), zkmax)
          cnew = c * exp((znew-z)*ff)
          !
          z    = znew
          c    = max(1.0e-6_fp, cnew)
       enddo
       aks_ss3d  = zkmax
       caks_ss3d = c
    endif
end subroutine calseddf2004

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
module bgrid_cold_start_mod

use    bgrid_horiz_mod, only: horiz_grid_type, get_horiz_grid_bound, TGRID
use     bgrid_halo_mod, only: update_halo, TEMP
use gaussian_topog_mod, only: gaussian_topog_init
use            fms_mod, only: file_exist, open_namelist_file, open_restart_file, &
                              open_ieee32_file, check_nml_error, close_file,     &
                              mpp_pe, mpp_root_pe, set_domain, read_data,        &
                              stdlog, error_mesg, FATAL, write_version_number,   &
                              field_size, uppercase
use      constants_mod, only: GRAV, RDGAS

implicit none
private

public cold_start_resol, cold_start

!-----------------------------------------------------------------------
!----- namelist ------
!
!  nlon = number of grid points along the longitude axis (1st dimension)
!  nlat = number of grid points along the latitude axis (2nd dimension)
!  nlev = number of vertical levels (equally spaced in sigma)
!  pref = initial surface pressure in pascals
!  tref = initial temperature in deg kelvin
!  equal_vert_spacing = specifies whether equal vertical spacing in sigma
!                       should be used (TRUE) or unequal spacing using
!                       the Smagorinski formula (FALSE).
!  topog_option = how to compute topography
!                  possible values are: FLAT, FILE, GAUSS
!  topog_file   = name of topography restart file in topog_option=FILE
!
!  NOTE: nlon and nlat are for the global compute grid
!

integer :: nlon = 0
integer :: nlat = 0
integer :: nlev = 0
   real :: pref = 1000.e2
   real :: tref = 255.
logical :: equal_vert_spacing = .true.
character(len=8)   :: topog_option = 'FLAT'
character(len=128) :: topog_file = ' '

namelist /bgrid_cold_start_nml/ nlon, nlat, nlev, pref, tref, &
                                equal_vert_spacing, topog_option, &
                                topog_file

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

character(len=128) :: version = '$Id: bgrid_cold_start.f90,v 11.0 2004/09/28 19:07:13 fms Exp $'
character(len=128) :: tag = '$Name: latest $'

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

contains

!#######################################################################

subroutine cold_start_resol ( nx, ny, nz )

      integer, intent(out) :: nx, ny, nz

      integer :: unit, ierr, io

!-------------- read namelist --------------

   if ( file_exist('input.nml')) then
      unit = open_namelist_file ( )
      ierr=1; do while (ierr /= 0)
         read  (unit, nml=bgrid_cold_start_nml, iostat=io, end=10)
         ierr = check_nml_error(io,'bgrid_cold_start_nml')
      enddo
 10   call close_file (unit)
   endif

!-------- write version and namelist to log file --------

   call write_version_number (version,tag)
   if (mpp_pe() == mpp_root_pe()) write (stdlog(), nml=bgrid_cold_start_nml)

!------- must specify a resolution -----

   if (nlon == 0 .or. nlat == 0 .or. nlev == 0)  &
   call error_mesg ('bgrid_cold_start_mod', 'resolution not specified', FATAL)

!------- otherwise, return resolution to calling program ------

   nx  = nlon
   ny  = nlat
   nz  = nlev

end subroutine cold_start_resol

!#######################################################################

subroutine cold_start ( Hgrid, eta, peta, fis, res, ps, pssl, u, v, t )

   type(horiz_grid_type), intent(inout) :: Hgrid
   real, intent(out), dimension(:)      :: eta, peta
   real, intent(out), dimension(Hgrid%ilb:,Hgrid%jlb:)     :: fis, res, &
                                                              ps, pssl
   real, intent(out), dimension(Hgrid%ilb:,Hgrid%jlb:,:)   :: u, v, t
  !real, intent(out), dimension(Hgrid%ilb:,Hgrid%jlb:,:,:) :: r
   integer :: m
   real    :: rt

!   very simple initial condition
!   -----------------------------
!   uniform pressure (adjusted for topog)
!   isothermal
!   no wind

!--- no hybrid option ---
   peta = 0.0

   eta = compute_sigma (size(eta(:))-1)

   call compute_topog (Hgrid, fis, res)

   rt = RDGAS * tref

     ps   = pref*exp(-fis/rt)
     pssl = ps                ! no eta option
     u    = 0.0
     v    = 0.0
     t    = tref

end subroutine cold_start

!#######################################################################

 subroutine compute_topog (Hgrid, fis, res)

   type(horiz_grid_type), intent(inout) :: Hgrid
   real, intent(out), dimension(Hgrid%ilb:,Hgrid%jlb:) :: fis, res

   integer :: unit, ires, jres, siz(4), nc
   logical :: eof_found
!-----------------------------------------------------------------------

!---------------------------------------------
! OPTION 1: No mountains, flat topography
!---------------------------------------------
  if ( trim(uppercase(topog_option)) == 'FLAT' ) then
       fis = 0.
       res = 1.

!---------------------------------------------
! OPTION 2: gaussian-shaped mountains
!---------------------------------------------
  else if ( trim(uppercase(topog_option(1:5))) == 'GAUSS' ) then

       call gaussian_topog_init ( Hgrid%Tmp%alm(:,Hgrid%Tmp%js), &
                                  Hgrid%Tmp%aph(Hgrid%Tmp%is,:), &
                                  fis(:,:) )
       fis = fis * GRAV                       ! convert topog from meters to geop ht
       where (fis(:,:) < 0.0) fis(:,:) = 0.0  ! only in case of roundoff error
       res = 1.

!---------------------------------------------
! OPTION 3: read topography from restart file
!---------------------------------------------
  else if ( trim(uppercase(topog_option)) == 'FILE' ) then

   if ( file_exist(trim(topog_file)) ) then
      ! netcdf or native?
        nc = len_trim(topog_file)
      !------------------------------------------
      ! read topography from netcdf restart file

        if (topog_file(nc-2:nc) == '.nc') then

          ! read fis
            call field_size (topog_file(1:nc), 'fis', siz )
            if ( siz(1) /= Hgrid%nlon .or. siz(2) /= Hgrid%nlat )     &
                 call error_mesg ('bgrid_cold_start_mod', 'incorrect resolution &
                              &or no fis field in netcdf topography file', FATAL)
            call read_data (topog_file(1:nc), 'fis', fis, Hgrid%Tmp%Domain)
            ! read res (if it exists)
            call field_size (topog_file(1:nc), 'res', siz )
            if ( siz(1) == Hgrid%nlon .and. siz(2) == Hgrid%nlat ) then
                call read_data (topog_file(1:nc), 'res', res, Hgrid%Tmp%Domain)
            else
                res = 1.
            endif
      !------------------------------------------
      ! read topography from native restart file

        else

            unit = open_restart_file (topog_file(1:nc), action='read')
            read (unit) ires,jres
            if ( ires /= Hgrid%nlon .or. jres /= Hgrid%nlat )  &
                            call error_mesg ('bgrid_cold_start_mod',  &
                        'incorrect resolution in native topography file', FATAL)
            call set_domain (Hgrid%Tmp%Domain)
            call read_data  (unit, fis)                ! read topog
            call read_data  (unit, res, end=eof_found) ! read res (if it exists)
            if (eof_found) res = 1.
            call close_file (unit)

        endif

   else
      ! restart file does not exist
        call error_mesg ('bgrid_cold_start',  &
                         'restart file for topography does not exist', FATAL)
   endif

!---------------------------------------------
! OPTION 4: error condition
!---------------------------------------------
  else
        call error_mesg ('bgrid_cold_start',  &
                         'invalid value for topog_option', FATAL)
  endif

! halos
  call update_halo ( Hgrid, TEMP, fis )
  call update_halo ( Hgrid, TEMP, res )

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

 end subroutine compute_topog

!#######################################################################

 function compute_sigma (nlev) result (eta)

   integer, intent(in) :: nlev
   real, dimension(nlev+1) :: eta
   real :: dz, qk
   integer :: k

     eta(1) = 0.0
     eta(nlev+1) = 1.0

   if ( equal_vert_spacing ) then

     ! levels with equal sigma spacing
       dz = 1./nlev
       do k = 2, nlev
         eta(k) = eta(k-1) + dz
       enddo

   else

     ! space levels according to Smagorinski (1965, MWR, pp.727-768)
     ! sigma spacing smaller at the top and bottom of model
       do k = 1, nlev-1
         qk = real(2*k)/real(2*nlev)
         eta(k+1) = qk*qk*(3.-2.*qk)
       enddo

   endif

 end function compute_sigma

!#######################################################################

end module bgrid_cold_start_mod

max_stars_repo_name: JoyMonteiro/FMSDynamicalCore


text:
CS    REAL FUNCTION EXPEI(X)
      DOUBLE PRECISION FUNCTION EXPEI(X)
C--------------------------------------------------------------------
C
C This function program computes approximate values for the
C   function  exp(-x) * Ei(x), where  Ei(x)  is the exponential
C   integral, and  x  is real.
C
C  Author: W. J. Cody
C
C  Latest modification: March 9, 1992
C
C--------------------------------------------------------------------
      INTEGER INT
CS    REAL  X, RESULT
      DOUBLE PRECISION  X, RESULT
C--------------------------------------------------------------------
      INT = 3
      CALL CALCEI(X,RESULT,INT)
      EXPEI = RESULT
      RETURN
C---------- Last line of EXPEI ----------
      END

ext: f

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
      program ecg
      call printout
      end
      subroutine printout
      print *, 'hello'
      end


text:
*+ iocmd_tclsetl

       subroutine iocmd_tclsetl( interp, val, var, el, s )
C      ---------------------------------------------------
C
C Set a logical value into a tcl variable
C
C Given:
C   pointer to interpreter
       integer    interp(*)
C   value to set
       logical    val
C   variable name
       character  var*(*)
C   element name if var is an array
       character  el*(*)
C
C Updated:
C   error status
       integer    s
C
C-
       include '/mrao/include/chrlib_functions.inc'
       character   name*64, str*256

       if (val) then
         str = '1'//char(0)
       else
         str = '0'//char(0)
       endif
       if (chr_lenb(el).gt.0) then
         name = var(1:chr_lenb(var))//
     *          '('//el(1:chr_lenb(el))//')'//char(0)
       else
         name = var(1:chr_lenb(var))//char(0)
       endif
       call iocmd_tclset( interp, name, str )

       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
Mason is currently the Administrative ViceChair of the California College Republicans (CCR). Officers of the Statewide College Republican federation are constitutionally exofficio officers in their local chapter club (Davis College Republicans).  Therefore Mason serves as a DCR and CCR officer simultaneously. Which also means
..Mason is a very busy man.

Most recently Mason was the Convention Director for the 2005 CCR Convention in Sacramento; this event gained attention on syndicated national radio programs such as Al Rantel and Mark Williams. The event is said to be the largest Republican youth gathering at the Capitol in State History, certainly a force to be reckoned with. He is best known on campus for videotaping the Gender and Sexuality Commissions Generation Sex Week. Although attorneys for Mason vouch for the legal soundness of this action, members who were upset by this claim that Mason violated safe zone regulations. Legal analysts continue to remind GASC members that there is no such thing as a safe zone described in the law, and as long as the event was paid for by student funds and opened to the public, it was entirely within Masons right to videotape the event.

Political Leadership Positions Held on Campus

  Administrative ViceChairman, California College Republicans (2005Present)

  Convention Director, California College Republicans (2005)

  Capitol Region ViceChairman, California College Republicans (200405)

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c***************************************************************************
c			OBL_ACC.F
c*************************************************************************
c OBL_ACC returns the BARYCENTRIC x,y,z components of the acc. on NBOD
c particles due to the oblateness of mass(1) using
c the values of J2RP2 and J4RP4 passed into the routine.
c (J2RP2 for example is the product of
c J_2 times the square of the central body's radius)
c Here we return the net acc. produced
c only by the J2 and J4 terms (i.e. including
c neither the monopole nor higher order terms).
c
c
c             Input:
c                 nbod     ==>  number of massive bodies (incl. central one)
c                 mass(*)  ==>  masses of particles (real*8 array)
c                 j2rp2    ==>  scaled value of j2 moment (real*8 scalar)
c                 j4rp4    ==>  scaled value of j4 moment (real*8 scalar)
c                                    (real*8 vectors)
c                 xh(*),yh(*),zh(*)   ==>  HELIO. positions of particles
c                 irh(*)   ==> 1./ magnitude of radius vector (real*8 vector)
c                                (passed in to save calcs.)
c             Output:
c               aoblx(*),aobly(*),aoblz(*)  ==>  BARY. components of accel
c                                        (real*8 vectors)
c
c Remarks:  aoblx(1) (for example) contains x-component of
c           bary. acc. of central body
c Authors:  Martin Duncan
c Date:    3/4/94
c Last revision:

      subroutine obl_acc(nbod,mass,j2rp2,j4rp4,xh,yh,zh,irh,
     &     aoblx,aobly,aoblz)

      include '../swift.inc'

c...  Inputs Only:
      integer nbod
      real*8 j2rp2,j4rp4
      real*8 mass(NPLMAX)
      real*8 xh(NPLMAX),yh(NPLMAX),zh(NPLMAX),irh(NPLMAX)

c...  Output
      real*8 aoblx(NPLMAX),aobly(NPLMAX),aoblz(NPLMAX)

c...  Internals
      integer n
      real*8 rinv2,t0,t1,t2,t3
      real*8 fac1,fac2

c----
c...  executable code

c First get the bary acc. of each "planet" due to central oblate "sun"

	do n =2,nbod

c Note that here we assume we know inverse of radius rather than calc. it
c from (x,y,z) to save the sqrt.
	  rinv2= irh(n)**2
	  t0 = -mass(1)*rinv2*rinv2*irh(n)
	  t1 = 1.5d0 *j2rp2
	  t2 = zh(n)*zh(n)*rinv2
	  t3 = 1.875d0 *j4rp4*rinv2

	  fac1 = t0*(t1 - t3 - (5.d0*t1 - (14.d0 - 21.d0*t2)*t3)*t2)
	  fac2 = 2.d0*t0*(t1 - (2.d0 - (14.d0*t2/3.d0))*t3)

          aoblx(n) = fac1*xh(n)
          aobly(n) = fac1*yh(n)
          aoblz(n) = (fac1 + fac2)*zh(n)

	enddo

c Now compute the bary. acc. of Sun due to all the planets

	aoblx(1) = 0.d0
	aobly(1) = 0.d0
	aoblz(1) = 0.d0
	do n=2,nbod
	   aoblx(1) = aoblx(1) - mass(n)*aoblx(n)/mass(1)
	   aobly(1) = aobly(1) - mass(n)*aobly(n)/mass(1)
	   aoblz(1) = aoblz(1) - mass(n)*aoblz(n)/mass(1)
	enddo

        return
        end                       !  obl_acc.f
c____________________________________________________________________________

max_stars_repo_name: Simske/exostriker


text:
module Mod_nsm_Viscosity
   use typre
   use Mod_Element
contains

   subroutine nsi_vislaw(e,grvel,lawvi,LawViParam,acvis)
      implicit none
      class(FiniteElement)     :: e
      integer(ip), intent(in)  :: lawvi
      real(rp),    intent(in)  :: LawViParam(*)
      real(rp),    intent(in)  :: grvel(e%ndime,e%ndime)
      real(rp),    intent(out) :: acvis
      integer(ip) :: idime,jdime,inode,vislaw
      real(rp)    :: vis_a,vispolaw,polawn,viszer,visinf,vislamb
      real(rp)    :: I2,strain(e%ndime,e%ndime),facvis1,facvis2,heavi

      !********************************************************************
      !vislawParam(1)=newtonian viscosity (mu) or consistency index (m)
      !visco(2)=infinity viscosity Carreau-Yasuda model
      !visco(3)=power law index in non newtonian models
      !visco(4)=lambda parameter in Carreau Yasuda model
      !visco(5)=a parameter in carreau yasuda model if a=2 carreau model
      !********************************************************************

      vispolaw=LawViParam(1)
      polawn=LawViParam(2)
      viszer=LawViParam(1)
      visinf=LawViParam(3)
      vislamb=LawViParam(4)
      vis_a=LawViParam(5)
      vislaw=lawvi

      !rate of strain tensor to use in non Newtonian models
      I2 = 0.0_rp
      strain = grvel + transpose(grvel)
      call tracedoubledot(e%ndime,strain,I2)
      !call trace(e%ndime,strain,tr_strain1)
      !call trace(e%ndime,matmul(strain,strain),tr_strain2)
      !I2 = tr_strain1**2 + tr_strain2

      if (vislaw == 1) then
         acvis = vispolaw*((0.5_rp*I2)**(polawn*0.5_rp-0.5_rp)) !Power Law viscosity
         if (polawn .gt. 1.0_rp) then
            facvis1 = LawViParam(1)/10000_rp
            facvis2 = LawViParam(1)*10000_rp
         end if
         if (polawn .le. 1.0_rp) then
            facvis1 = LawViParam(1)/10000_rp
            facvis2 = LawViParam(1)*10000_rp
         end if
         if (acvis <= facvis1) acvis = facvis1
         if (acvis >= facvis2) acvis = facvis2
      elseif (vislaw == 2) then
         acvis = visinf+(viszer-visinf)*(1.0_rp+(0.5_rp*vislamb*I2)**vis_a)**((polawn-1.0_rp)/vis_a)
      end if

   end subroutine nsi_vislaw

   subroutine nsm_smago(e,grvel,denac,turbu,vista)
      use typre
      use Mod_Element
      implicit none
      class(FiniteElement)  :: e
      real(rp), intent(in)  :: turbu
      real(rp), intent(in)  :: grvel(e%ndime,e%ndime)
      real(rp), intent(in)  :: denac
      real(rp), intent(out) :: vista
      real(rp)    :: seci4,divel
      integer(ip) :: igaus,inode,idime,jdime

      vista = 0.0_rp
      !Smagorinsky
      seci4=0.0_rp
      do idime=1,e%ndime
         do jdime=1,e%ndime
            seci4=seci4+grvel(idime,jdime)*(grvel(idime,jdime)+grvel(jdime,idime))
         end do
      end do
      divel=0.0_rp
      do idime=1,e%ndime
         divel=divel+grvel(idime,idime)
      end do
      vista=2.0_rp*denac*turbu*e%hleng(e%ndime)*e%hleng(e%ndime)*sqrt(seci4-2.0_rp*divel*divel/3.0_rp)

   end subroutine nsm_smago

   !> This is the WALE LES model based on  Subgrid-Scale Stress Modelling Based on the
   !>   Square of the Velocity Gradient Tensor. F. NICOUD and F. DUCROS Flow, Turbulence and Combustion 62: 183–200, 1999.
   subroutine nsm_wale(e,grvel,denac,turbu,vista)
      use typre
      use Mod_Element
      implicit none
      class(FiniteElement)  :: e
      real(rp), intent(in)  :: turbu
      real(rp), intent(in)  :: grvel(e%ndime,e%ndime)
      real(rp), intent(in)  :: denac
      real(rp), intent(out) :: vista
      real(rp) :: CaligraphicS(e%ndime,e%ndime)
      real(rp) :: S(e%ndime,e%ndime), g2(e%ndime,e%ndime)
      real(rp) :: MeanStressG2, OP1,OP2, TraceG2
      real(rp) :: CSDoubleContractionCS, SDoubleContractionS
      integer(ip) :: idime,jdime

      !Default value
      vista = 0.0_rp
      !Symmetric gradient
      S = 0.5_rp*(grvel + transpose(grvel))
      !Gradient times gradient
      g2 = matmul(grvel,grvel)

      TraceG2 = 0.0_rp
      do idime = 1,e%ndime
         TraceG2 = TraceG2 + g2(idime,idime)
      enddo
      MeanStressG2 = 1.0_rp/3.0_rp*TraceG2

      !CaligraphicS
      CaligraphicS = 0.5_rp*(g2 + transpose(g2))
      do idime = 1,e%ndime
         CaligraphicS(idime,idime) = CaligraphicS(idime,idime) - MeanStressG2
      enddo

      !Symmetric gradient against itself
      SDoubleContractionS = 0.0_rp
      do idime = 1,e%ndime
         do jdime = 1,e%ndime
            SDoubleContractionS = SDoubleContractionS + S(idime,jdime)*S(idime,jdime)
         enddo
      enddo

      !Caligraphic S against itself
      CSDoubleContractionCS = 0.0_rp
      do idime = 1,e%ndime
         do jdime = 1,e%ndime
            CSDoubleContractionCS = CSDoubleContractionCS + CaligraphicS(idime,jdime)*CaligraphicS(idime,jdime)
         enddo
      enddo

      !OP1 and OP2 as defined in the Nicoud paper
      OP1 = CSDoubleContractionCS**(3.0_rp/2.0_rp)
      OP2 = SDoubleContractionS**(5.0_rp/2.0_rp) + CSDoubleContractionCS**(5.0_rp/4.0_rp)

      if (OP2 /= 0.0_rp) vista = denac*((turbu*e%hleng(e%ndime))**2.0_rp)*OP1/OP2

   end subroutine nsm_wale

end module

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM