New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
usrdef_istate.F90 in NEMO/branches/UKMO/NEMO_4.0.4_hpge_ovf/src/OCE/USR – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_hpge_ovf/src/OCE/USR/usrdef_istate.F90 @ 15676

Last change on this file since 15676 was 15676, checked in by dbruciaferri, 2 years ago

Adding Shchepetkin & McWilliams? 2003 initial profile
File size: 10.9 KB
Line 
1MODULE usrdef_istate
2   !!======================================================================
3   !!                   ***  MODULE  usrdef_istate   ***
4   !!
5   !!                     ===  GYRE configuration  ===
6   !!
7   !! User defined : set the initial state of a user configuration
8   !!======================================================================
9   !! History :  4.0 ! 2016-03  (S. Flavoni) Original code
10   !!                ! 2021-10  (D. Bruciaferri) JMMP initial condition for
11   !!                                            HPGE test
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!  usr_def_istate : initial state in Temperature and salinity
16   !!----------------------------------------------------------------------
17   USE par_oce        ! ocean space and time domain
18   USE phycst         ! physical constants
19   USE splines
20   USE dtatsd         ! data temperature and salinity   (dta_tsd routine)
21   !
22   USE in_out_manager ! I/O manager
23   USE lib_mpp        ! MPP library
24   
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   usr_def_istate   ! called in istate.F90
29
30   !!----------------------------------------------------------------------
31   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
32   !! $Id: usrdef_istate.F90 10069 2018-08-28 14:12:24Z nicolasmartin $
33   !! Software governed by the CeCILL license (see ./LICENSE)
34   !!----------------------------------------------------------------------
35CONTAINS
36 
37   SUBROUTINE usr_def_istate( pdept, ptmask, pts, pu, pv, pssh )
38      !!----------------------------------------------------------------------
39      !!                   ***  ROUTINE usr_def_istate  ***
40      !!
41      !! ** Purpose :   Initialization of the dynamics and tracers
42      !!                Here GYRE configuration example : (double gyre with rotated domain)
43      !!
44      !! ** Method  : - set temprature field
45      !!              - set salinity   field
46      !!----------------------------------------------------------------------
47      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pdept   ! depth of t-point               [m]
48      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   ptmask  ! t-point ocean mask             [m]
49      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pts     ! T & S fields      [Celsius ; g/kg]
50      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pu      ! i-component of the velocity  [m/s]
51      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pv      ! j-component of the velocity  [m/s]
52      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(  out) ::   pssh    ! sea-surface height
53      !
54      REAL(wp), DIMENSION(42)  ::   zdep, zsal, ztem
55      REAL(wp) :: a0, a1, a2, b0, b1, b2
56      INTEGER :: ji, jj, jk  ! dummy loop indices
57      !!----------------------------------------------------------------------
58      !
59      IF(lwp) WRITE(numout,*)
60      IF(lwp) WRITE(numout,*) 'usr_def_istate : analytical definition of initial state '
61      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   Ocean at rest, with an horizontally uniform T and S profiles'
62      !
63      pu  (:,:,:) = 0._wp        ! ocean at rest
64      pv  (:,:,:) = 0._wp
65      pssh(:,:)   = 0._wp
66      !
67     
68      zdep(:) = (/ &
69                   &    5.02,   15.08,   25.16,   35.28,   45.45,   55.69, &
70                   &   66.04,   76.55,   87.27,   98.31,  109.81,  121.95, &
71                   &  135.03,  149.43,  165.73,  184.70,  207.43,  235.39, &
72                   &  270.53,  315.37,  372.97,  446.80,  540.50,  657.32, &
73                   &  799.55,  968.00, 1161.81, 1378.66, 1615.29, 1868.07, &
74                   & 2133.52, 2408.58, 2690.78, 2978.17, 3269.28, 3563.04, &
75                   & 3858.68, 4155.63, 4453.50, 4752.02, 5050.99, 6300.27  &
76                   & /)
77
78      SELECT CASE(nn_tsd_type)
79     
80         CASE(1) ! T/S July on-shelf
81            ztem(:) = (/ &
82                         & 13.8059, 13.1661, 12.0471, 10.7265,  9.5585,  8.7273, &
83                         &  8.2527,  7.9696,  7.7605,  7.5965,  7.5038,  7.5018, &
84                         &  7.5609,  7.6436,  7.7050,  7.7011,  7.5887,  7.3471, &
85                         &  6.9203,  6.1857,  5.2453,  4.1932,  3.4076,  2.8450, &
86                         &  2.6111,  2.5000,  2.5000,  2.5000,  2.5000,  2.5000, &
87                         &  2.5000,  2.5000,  2.5000,  2.5000,  2.5000,  2.5000, &
88                         &  2.5000,  2.5000,  2.5000,  2.5000,  2.5000,  2.5000  &
89                         & /)
90
91            zsal(:) = (/ &
92                         & 33.7727, 33.9841, 34.2683, 34.5326, 34.7335, 34.9079, &
93                         & 35.0119, 35.0726, 35.0902, 35.1034, 35.1101, 35.1223, &
94                         & 35.1436, 35.1718, 35.1893, 35.1929, 35.1828, 35.1670, &
95                         & 35.1408, 35.0990, 35.0415, 34.9760, 34.9168, 34.8753, &
96                         & 34.8558, 34.8502, 34.8504, 34.8506, 34.8508, 34.8510, &
97                         & 34.8512, 34.8514, 34.8516, 34.8518, 34.8520, 34.8522, &
98                         & 34.8524, 34.8526, 34.8528, 34.8530, 34.8532, 34.8534  &
99                         & /)
100
101         CASE(2) ! T/S July off-shelf
102     
103            ztem(:) = (/ &
104                         & 13.0669, 12.8587, 12.4760, 11.9986, 11.5363, 11.1627, &
105                         & 10.8898, 10.6753, 10.4927, 10.3334, 10.2182, 10.1457, &
106                         & 10.1038, 10.0734, 10.0389,  9.9968,  9.9459,  9.8836, &
107                         &  9.8069,  9.6953,  9.5345,  9.2901,  8.9319,  8.4192, &
108                         &  7.7006,  6.7895,  5.7774,  4.8576,  4.1510,  3.6716, &
109                         &  3.3331,  3.0606,  2.8275,  2.6317,  2.4735,  2.3497, &
110                         &  2.2601,  2.1973,  2.1555,  2.1237,  2.1072,  2.1000  &
111                         & /)
112
113            zsal(:) = (/ &
114                         & 35.2001, 35.2052, 35.2186, 35.2411, 35.2661, 35.2873, &
115                         & 35.3021, 35.3124, 35.3205, 35.3267, 35.3304, 35.3330, &
116                         & 35.3355, 35.3393, 35.3422, 35.3438, 35.3436, 35.3428, &
117                         & 35.3413, 35.3374, 35.3313, 35.3239, 35.3192, 35.3171, &
118                         & 35.3171, 35.3171, 35.3171, 35.3171, 35.3171, 35.3171, &
119                         & 35.3171, 35.3171, 35.3171, 35.3171, 35.3171, 35.3171, &
120                         & 35.3171, 35.3171, 35.3171, 35.3171, 35.3171, 35.3171  &
121                         & /)
122
123         CASE(3) ! T/S January on-shelf
124
125            ztem(:) = (/ &
126                         &  7.8383,  7.9482,  8.0837,  8.1641,  8.1750,  8.1510, &
127                         &  8.1110,  8.0510,  7.9648,  7.8804,  7.8221,  7.8020, &
128                         &  7.7959,  7.7919,  7.8288,  7.8825,  7.9012,  7.7028, &
129                         &  7.2095,  6.3642,  5.2466,  4.0221,  2.9315,  2.1751, &
130                         &  1.7760,  1.6325,  1.6008,  1.6000,  1.6000,  1.6000, &
131                         &  1.6000,  1.6000,  1.6000,  1.6000,  1.6000,  1.6000, &
132                         &  1.6000,  1.6000,  1.6000,  1.6000,  1.6000,  1.6000  &
133                         & /)
134
135            zsal(:) = (/ &
136                         & 34.0842, 34.1682, 34.3100, 34.4740, 34.6210, 34.7275, &
137                         & 34.7861, 34.8231, 34.8661, 34.9293, 34.9936, 35.0405, &
138                         & 35.0618, 35.0683, 35.0791, 35.0980, 35.1224, 35.1287, &
139                         & 35.1140, 35.0773, 35.0315, 34.9838, 34.9431, 34.9159, &
140                         & 34.9019, 34.8971, 34.8960, 34.8962, 34.8964, 34.8966, &
141                         & 34.8968, 34.8970, 34.8972, 34.8974, 34.8976, 34.8978, &
142                         & 34.8980, 34.8982, 34.8984, 34.8986, 34.8988, 34.8990  &
143                         & /)
144
145         CASE(4) ! T/S January off-shelf
146
147            ztem(:) = (/ &
148                         &  9.9555,  9.9499,  9.9411,  9.9311,  9.9211,  9.9111, &
149                         &  9.9011,  9.8911,  9.8811,  9.8711,  9.8611,  9.8511, &
150                         &  9.8411,  9.8311,  9.8211,  9.8111,  9.8011,  9.7911, &
151                         &  9.7840,  9.7487,  9.6420,  9.3992,  9.0030,  8.4382, &
152                         &  7.7101,  6.8430,  5.9185,  5.0553,  4.3437,  3.8076, &
153                         &  3.4096,  3.0998,  2.8462,  2.6397,  2.4862,  2.3868, &
154                         &  2.3374,  2.3159,  2.3059,  2.2959,  2.2870,  2.2815  &
155                         & /)
156
157            zsal(:) = (/ &
158                         & 35.2802, 35.2808, 35.2817, 35.2827, 35.2837, 35.2847, &
159                         & 35.2857, 35.2867, 35.2877, 35.2887, 35.2897, 35.2907, &
160                         & 35.2917, 35.2927, 35.2937, 35.2947, 35.2957, 35.2967, &
161                         & 35.2998, 35.3032, 35.3048, 35.2985, 35.2851, 35.2640, &
162                         & 35.2340, 35.1899, 35.1320, 35.0692, 35.0150, 34.9784, &
163                         & 34.9598, 34.9527, 34.9502, 34.9474, 34.9436, 34.9397, &
164                         & 34.9379, 34.9378, 34.9388, 34.9398, 34.9407, 34.9413  &
165                         & /)
166         CASE(5) ! Analytic T/S July on-shelf
167           ! a0 = 13.0
168           ! a1 = 0.00507078656
169           ! a2 = 2.37539619
170           ! b0 = -1.1
171           ! b1 = 0.01
172           ! b2 = 34.85
173         CASE(6) ! Analytic T in the tropical Atlantic ocean
174            rn_sal_sf = 35
175            a0 = 0.2 
176            a1 = 6.0 
177            a2 = 20. 
178            b0 = 2500.
179            b1 = 250.
180      END SELECT
181
182      IF (nn_tsd_type < 5) THEN  ! Use spline smoothing
183                                 ! horizontally uniform T & S profiles
184         DO jj = 1, jpj
185            DO ji = 1, jpi
186               pts(ji,jj,:,jp_tem) = spline3(zdep,ztem,pdept(ji,jj,:)) !* ptmask(ji,jj,:)
187               pts(ji,jj,:,jp_sal) = spline3(zdep,zsal,pdept(ji,jj,:)) !* ptmask(ji,jj,:)
188            END DO
189         END DO
190
191      ELSE IF (nn_tsd_type == 5) THEN ! Siddorn & Furner 2013 analytical function
192
193         DO jk = 1, jpk
194            DO jj = 1, jpj
195               DO ji = 1, jpi
196                  pts(ji,jj,jk,jp_sal) = rn_sal_sf !* ptmask(ji,jj,jk)
197                  pts(ji,jj,jk,jp_tem) = ( rn_c0_sf * ( 1-TANH(((pdept(ji,jj,jk)-rn_mld_sf)/20)*3.1415927/180) ) &
198                                         + rn_c1_sf * ((rn_maxdep_sf- pdept(ji,jj,jk)) / rn_maxdep_sf) ) !* ptmask(ji,jj,jk)
199               END DO
200            END DO
201         END DO
202
203      ELSE IF (nn_tsd_type == 6) THEN ! Shchepetkin and McWilliams 2003 analytical function
204
205         DO jk = 1, jpk
206            DO jj = 1, jpj
207               DO ji = 1, jpi
208                  pts(ji,jj,jk,jp_sal) = rn_sal_sf !* ptmask(ji,jj,jk)
209                  pts(ji,jj,jk,jp_tem) = a0 + a1*EXP(-pdept(ji,jj,jk)/b0) &
210                                            + a2*EXP(-pdept(ji,jj,jk)/b1) !* ptmask(ji,jj,jk)
211               END DO
212            END DO
213         END DO
214
215      END IF
216
217   END SUBROUTINE usr_def_istate
218
219   !!======================================================================
220END MODULE usrdef_istate
Note: See TracBrowser for help on using the repository browser.