[6331] | 1 | MODULE function_ana |
---|
| 2 | ! |
---|
| 3 | IMPLICIT NONE |
---|
| 4 | ! |
---|
| 5 | CONTAINS |
---|
| 6 | ! |
---|
| 7 | !**************************************************************************************** |
---|
| 8 | SUBROUTINE function_ana1(ni, nj, xcoor, ycoor, fnc_ana) |
---|
| 9 | !**************************************************************************************** |
---|
| 10 | |
---|
| 11 | !**** *function ana* - calculate analytical function |
---|
| 12 | |
---|
| 13 | INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) ! double |
---|
| 14 | |
---|
| 15 | INTEGER, INTENT(IN) :: ni, nj |
---|
| 16 | REAL(kind=wp), DIMENSION(ni,nj), INTENT(IN) :: xcoor, ycoor |
---|
| 17 | REAL(kind=wp), DIMENSION(ni,nj), INTENT(OUT) :: fnc_ana |
---|
| 18 | |
---|
| 19 | REAL (kind=wp), PARAMETER :: dp_pi=3.14159265359 |
---|
| 20 | REAL (kind=wp), PARAMETER :: dp_conv = dp_pi/180. |
---|
| 21 | REAL(kind=wp) :: dp_length, coef, coefmult |
---|
| 22 | INTEGER :: i,j |
---|
| 23 | CHARACTER(LEN=7) :: cl_anaftype="fcos" |
---|
| 24 | |
---|
| 25 | DO j=1,nj |
---|
| 26 | DO i=1,ni |
---|
| 27 | |
---|
| 28 | SELECT CASE (cl_anaftype) |
---|
| 29 | CASE ("fcos") |
---|
| 30 | dp_length = 1.2*dp_pi |
---|
| 31 | coef = 2. |
---|
| 32 | coefmult = 1. |
---|
| 33 | fnc_ana(i,j) = coefmult*(coef - COS( dp_pi*(ACOS( COS(xcoor(i,j)*dp_conv)*COS(ycoor(i,j)*dp_conv) )/dp_length)) ) |
---|
| 34 | |
---|
| 35 | CASE ("fcossin") |
---|
| 36 | dp_length = 1.d0*dp_pi |
---|
| 37 | coef = 21.d0 |
---|
| 38 | coefmult = 3.846d0 * 20.d0 |
---|
| 39 | fnc_ana(i,j) = coefmult*(coef - COS( dp_pi*(ACOS( COS(ycoor(i,j)*dp_conv)*COS(ycoor(i,j)*dp_conv) )/dp_length)) * & |
---|
| 40 | SIN( dp_pi*(ASIN( SIN(xcoor(i,j)*dp_conv)*SIN(ycoor(i,j)*dp_conv) )/dp_length)) ) |
---|
| 41 | END SELECT |
---|
| 42 | |
---|
| 43 | ENDDO |
---|
| 44 | ENDDO |
---|
| 45 | |
---|
| 46 | END SUBROUTINE function_ana1 |
---|
| 47 | |
---|
| 48 | !!**************************************************************************************** |
---|
| 49 | !SUBROUTINE function_ana2(lon, lat, mask, fnc_ana) |
---|
| 50 | SUBROUTINE function_ana2(ni, nj, lon, lat, fnc_ana) |
---|
| 51 | !!**************************************************************************************** |
---|
| 52 | |
---|
| 53 | !**** *function ana* - calculate analytical function |
---|
| 54 | |
---|
| 55 | INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) ! double |
---|
| 56 | |
---|
| 57 | INTEGER, INTENT(IN) :: ni, nj |
---|
| 58 | REAL (kind=wp), DIMENSION(ni,nj), INTENT(IN) :: lon, lat |
---|
| 59 | REAL(kind=wp), DIMENSION(ni,nj), INTENT(OUT) :: fnc_ana |
---|
| 60 | |
---|
| 61 | REAL (kind=wp), PARAMETER :: coef=2., dp_pi=3.14159265359 |
---|
| 62 | REAL (kind=wp) :: dp_length, dp_conv |
---|
| 63 | INTEGER :: i, j |
---|
| 64 | |
---|
| 65 | ! Analytical Gulf Stream |
---|
| 66 | REAL (kind=wp) :: gf_coef, gf_ori_lon, gf_ori_lat, & |
---|
| 67 | & gf_end_lon, gf_end_lat, gf_dmp_lon, gf_dmp_lat |
---|
| 68 | REAL (kind=wp) :: gf_per_lon |
---|
| 69 | REAL (kind=wp) :: dx, dy, dr, dth, dc, dr0, dr1 |
---|
| 70 | |
---|
| 71 | ! Parameters for analytical function |
---|
| 72 | dp_length = 1.2*dp_pi |
---|
| 73 | dp_conv = dp_pi/180. |
---|
| 74 | gf_coef = 1.0 ! Coefficient for Gulf Stream term (0.0 = no Gulf Stream) |
---|
| 75 | gf_ori_lon = -80.0 ! Origin of the Gulf Stream (longitude in deg) |
---|
| 76 | gf_ori_lat = 25.0 ! Origin of the Gulf Stream (latitude in deg) |
---|
| 77 | gf_end_lon = -1.8 ! End point of the Gulf Stream (longitude in deg) |
---|
| 78 | gf_end_lat = 50.0 ! End point of the Gulf Stream (latitude in deg) |
---|
| 79 | gf_dmp_lon = -25.5 ! Point of the Gulf Stream decrease (longitude in deg) |
---|
| 80 | gf_dmp_lat = 55.5 ! Point of the Gulf Stream decrease (latitude in deg) |
---|
| 81 | |
---|
| 82 | ! ALLOCATE(fnc_ana(0:xy-1)) |
---|
| 83 | |
---|
| 84 | dr0 = SQRT(((gf_end_lon - gf_ori_lon)*dp_conv)**2 + & |
---|
| 85 | & ((gf_end_lat - gf_ori_lat)*dp_conv)**2) |
---|
| 86 | dr1 = SQRT(((gf_dmp_lon - gf_ori_lon)*dp_conv)**2 + & |
---|
| 87 | & ((gf_dmp_lat - gf_ori_lat)*dp_conv)**2) |
---|
| 88 | |
---|
| 89 | DO j=1,nj |
---|
| 90 | DO i=1,ni |
---|
| 91 | |
---|
| 92 | ! Original OASIS fcos analytical test function |
---|
| 93 | fnc_ana(i,j)=(coef-COS(dp_pi*(ACOS(COS(lat(i,j)*dp_conv)*& |
---|
| 94 | & COS(lon(i,j)*dp_conv))/dp_length))) |
---|
| 95 | gf_per_lon = lon(i,j) |
---|
| 96 | IF (gf_per_lon > 180.0) gf_per_lon = gf_per_lon - 360.0 |
---|
| 97 | IF (gf_per_lon < -180.0) gf_per_lon = gf_per_lon + 360.0 |
---|
| 98 | dx = (gf_per_lon - gf_ori_lon)*dp_conv |
---|
| 99 | dy = (lat(i,j) - gf_ori_lat)*dp_conv |
---|
| 100 | dr = SQRT(dx*dx + dy*dy) |
---|
| 101 | dth = ATAN2(dy, dx) |
---|
| 102 | dc = 1.3*gf_coef |
---|
| 103 | IF (dr > dr0) dc = 0.0 |
---|
| 104 | IF (dr > dr1) dc = dc * COS(dp_pi*0.5*(dr-dr1)/(dr0-dr1)) |
---|
| 105 | fnc_ana(i,j) = fnc_ana(i,j) + & |
---|
| 106 | & (MAX(1000.0*SIN(0.4*(0.5*dr+dth) + & |
---|
| 107 | & 0.007*COS(50.0*dth) + 0.37*dp_pi),999.0) - 999.0) * dc |
---|
| 108 | |
---|
| 109 | ENDDO |
---|
| 110 | ENDDO |
---|
| 111 | |
---|
| 112 | END SUBROUTINE function_ana2 |
---|
| 113 | |
---|
| 114 | !**************************************************************************************** |
---|
| 115 | SUBROUTINE function_vortex(ni, nj, xcoor, ycoor, fnc_ana) |
---|
| 116 | !**************************************************************************************** |
---|
| 117 | !!! VORTEX FROM TEMPEST-REMAP (as in XIOS) |
---|
| 118 | ! |
---|
| 119 | INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) ! double |
---|
| 120 | ! |
---|
| 121 | INTEGER, INTENT(IN) :: ni, nj |
---|
| 122 | REAL(kind=wp), DIMENSION(ni,nj), INTENT(IN) :: xcoor, ycoor |
---|
| 123 | REAL(kind=wp), DIMENSION(ni,nj), INTENT(OUT) :: fnc_ana |
---|
| 124 | |
---|
| 125 | REAL(kind=wp), PARAMETER :: dp_pi = 3.14159265359 |
---|
| 126 | REAL(kind=wp), PARAMETER :: dLon0 = 5.5 |
---|
| 127 | REAL(kind=wp), PARAMETER :: dLat0 = 0.2 |
---|
| 128 | REAL(kind=wp), PARAMETER :: dR0 = 3.0 |
---|
| 129 | REAL(kind=wp), PARAMETER :: dD = 5.0 |
---|
| 130 | REAL(kind=wp), PARAMETER :: dT = 6.0 |
---|
| 131 | REAL(kind=wp) :: dp_length, dp_conv |
---|
| 132 | |
---|
| 133 | REAL(kind=wp) :: dSinC, dCosC, dCosT, dSinT |
---|
| 134 | REAL(kind=wp) :: dTrm, dX, dY, dZ |
---|
| 135 | REAL(kind=wp) :: dlon, dlat |
---|
| 136 | REAL(kind=wp) :: dRho, dVt, dOmega |
---|
| 137 | |
---|
| 138 | INTEGER :: i,j |
---|
| 139 | CHARACTER(LEN=7) :: cl_anaftype="vortex" |
---|
| 140 | |
---|
| 141 | dp_conv = dp_pi/180. |
---|
| 142 | dSinC = SIN( dLat0 ) |
---|
| 143 | dCosC = COS( dLat0 ) |
---|
| 144 | |
---|
| 145 | DO j=1,nj |
---|
| 146 | DO i=1,ni |
---|
| 147 | ! Find the rotated longitude and latitude of a point on a sphere |
---|
| 148 | ! with pole at (dLon0, dLat0). |
---|
| 149 | dCosT = COS( ycoor(i,j)*dp_conv ) |
---|
| 150 | dSinT = SIN( ycoor(i,j)*dp_conv ) |
---|
| 151 | |
---|
| 152 | dTrm = dCosT * COS( xcoor(i,j)*dp_conv - dLon0 ) |
---|
| 153 | dX = dSinC * dTrm - dCosC * dSinT |
---|
| 154 | dY = dCosT * SIN( xcoor(i,j)*dp_conv - dLon0 ) |
---|
| 155 | dZ = dSinC * dSinT + dCosC * dTrm |
---|
| 156 | |
---|
| 157 | dlon = ATAN2( dY, dX ) |
---|
| 158 | IF( dlon < 0.0 ) dlon = dlon + 2.0 * dp_pi |
---|
| 159 | dlat = ASIN( dZ ) |
---|
| 160 | |
---|
| 161 | dRho = dR0 * COS(dlat) |
---|
| 162 | dVt = 3.0 * SQRT(3.0)/2.0/COSH(dRho)/COSH(dRho)*TANH(dRHo) |
---|
| 163 | IF (dRho == 0.0) THEN |
---|
| 164 | dOmega = 0.0 |
---|
| 165 | ELSE |
---|
| 166 | dOmega = dVt / dRho |
---|
| 167 | END IF |
---|
| 168 | |
---|
| 169 | fnc_ana(i,j) = 2.0 * ( 1.0 + TANH( dRho / dD * SIN( dLon - dOmega * dT ) ) ) |
---|
| 170 | |
---|
| 171 | END DO |
---|
| 172 | END DO |
---|
| 173 | ! |
---|
| 174 | END SUBROUTINE function_vortex |
---|
| 175 | |
---|
| 176 | END MODULE function_ana |
---|