source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/MY_SRC/usrdef_istate.F90 @ 12443

Last change on this file since 12443 was 12443, checked in by davestorkey, 9 months ago

2020/KERNEL-03_Storkey_Coward_RK3_stage2: More variable renaming:
atfp → rn_atfp (use namelist parameter everywhere)
rdtbt → rDt_e
nn_baro → nn_e
rn_scal_load → rn_load
rau0 → rho0

  • Property svn:keywords set to Id
File size: 12.2 KB
Line 
1MODULE usrdef_istate
2   !!======================================================================
3   !!                     ***  MODULE usrdef_istate   ***
4   !!
5   !!                      ===  CANAL configuration  ===
6   !!
7   !! User defined : set the initial state of a user configuration
8   !!======================================================================
9   !! History :  NEMO ! 2017-11  (J. Chanut) Original code
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!  usr_def_istate : initial state in Temperature and salinity
14   !!----------------------------------------------------------------------
15   USE par_oce        ! ocean space and time domain
16   USE dom_oce       
17   USE phycst         ! physical constants
18   !
19   USE in_out_manager ! I/O manager
20   USE lib_mpp        ! MPP library
21   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
22   !   
23   USE usrdef_nam
24   
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   usr_def_istate   ! called by istate.F90
29
30   !!----------------------------------------------------------------------
31   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
32   !! $Id$
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 CANAL configuration
43      !!
44      !! ** Method  :   Set a gaussian anomaly of pressure and associated
45      !!                geostrophic velocities
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      INTEGER  :: ji, jj, jk, jl  ! dummy loop indices
55      REAL(wp) :: zx, zy, zP0, zumax, zlambda, zr_lambda2, zn2, zf0, zH, zrho1, za, zf, zdzF
56      REAL(wp) :: zpsurf, zdyPs, zdxPs
57      REAL(wp) :: zdt, zdu, zdv
58      REAL(wp) :: zjetx, zjety, zbeta
59      REAL(wp), DIMENSION(jpi,jpj)  ::   zrandom
60      !!----------------------------------------------------------------------
61      !
62      IF(lwp) WRITE(numout,*)
63      IF(lwp) WRITE(numout,*) 'usr_def_istate : CANAL configuration, analytical definition of initial state'
64      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   '
65      !
66      IF (ln_sshnoise) CALL RANDOM_NUMBER(zrandom)
67      zjetx = ABS(rn_ujetszx)/2.
68      zjety = ABS(rn_ujetszy)/2.
69      !
70      SELECT CASE(nn_initcase)
71      CASE(0)    ! rest
72         
73         ! sea level:
74         pssh(:,:) = 0.
75         ! temperature:
76         pts(:,:,:,jp_tem) = 10._wp
77         ! salinity: 
78         pts(:,:,:,jp_sal) = 35._wp
79         ! velocities:
80         pu(:,:,:) = 0.
81         pv(:,:,:) = 0.
82         
83      CASE(1)    ! geostrophic zonal jet from -zjety to +zjety
84
85         ! sea level:
86         SELECT CASE( nn_fcase )
87         CASE(0)    ! f = f0
88            ! sea level: ssh = - fuy / g
89            WHERE( ABS(gphit) <= zjety )
90               pssh(:,:) = - ff_t(:,:) * rn_uzonal * gphit(:,:) * 1.e3 / grav
91            ELSEWHERE
92               pssh(:,:) = - ff_t(:,:) * rn_uzonal * SIGN(zjety, gphit(:,:)) * 1.e3 / grav
93            END WHERE
94         CASE(1)    ! f = f0 + beta*y
95            ! sea level: ssh = - u / g * ( fy + 0.5 * beta * y^2 )
96            zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra
97            WHERE( ABS(gphit) <= zjety )
98               pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
99            ELSEWHERE
100               pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3   &
101                  &                             + 0.5 * zbeta * zjety * zjety * 1.e6 )
102            END WHERE
103         END SELECT
104         ! temperature:
105         pts(:,:,:,jp_tem) = 10._wp
106         ! salinity: 
107         pts(:,:,jpk,jp_sal) = 0.
108         DO jk=1, jpkm1
109            pts(:,:,jk,jp_sal) = gphit(:,:)
110         END DO
111         ! velocities:
112         pu(:,:,:) = 0.
113         DO jk=1, jpkm1
114            WHERE( ABS(gphit) <= zjety ) pu(:,:,jk) = rn_uzonal
115         END DO
116         pv(:,:,:) = 0.
117         !                 
118      CASE(2)    ! geostrophic zonal current shear
119     
120         ! sea level:
121         SELECT CASE( nn_fcase )
122         CASE(0)    ! f = f0
123            ! sea level: ssh = - fuy / g
124            WHERE( ABS(gphit) <= zjety )
125               pssh(:,:) = - ff_t(:,:) * rn_uzonal * ABS(gphit(:,:)) * 1.e3 / grav
126            ELSEWHERE
127               pssh(:,:) = - ff_t(:,:) * rn_uzonal * zjety * 1.e3 / grav
128            END WHERE
129         CASE(1)    ! f = f0 + beta*y
130            ! sea level: ssh = - u / g * ( fy + 0.5 * beta * y^2 )
131            zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra
132            WHERE( ABS(gphit) <= zjety )
133               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   &
134                  &        * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
135            ELSEWHERE
136               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   &
137                  &        * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3 + 0.5 * zbeta * zjety * zjety * 1.e6 )
138            END WHERE
139         END SELECT
140         ! temperature:
141         pts(:,:,:,jp_tem) = 10._wp
142         ! salinity: 
143         pts(:,:,:,jp_sal) = 2.
144         DO jk=1, jpkm1
145            WHERE( ABS(gphiv) <= zjety ) pts(:,:,jk,jp_sal) = 2. + SIGN(1.,gphiv(:,:))
146         END DO
147         ! velocities:
148         pu(:,:,:) = 0.
149         DO jk=1, jpkm1
150            WHERE( ABS(gphiv) <= zjety ) pu(:,:,jk) = SIGN(rn_uzonal,gphit(:,:))*SIGN(1.,rn_uzonal)
151            WHERE( ABS(gphiv) == 0.    ) pu(:,:,jk) = 0. 
152         END DO
153         pv(:,:,:) = 0.
154         !                 
155      CASE(3)    ! gaussian zonal currant
156
157         ! zonal current
158         DO jk=1, jpkm1
159            ! gphit and lambda are both in km
160            pu(:,:,jk) = rn_uzonal * EXP( - 0.5 * gphit(:,:)**2 / rn_lambda**2 )
161         END DO
162         
163         ! sea level:
164         pssh(:,1) = - ff_t(:,1) / grav * pu(:,1,1) * e2t(:,1)
165         DO jl=1, jpnj
166            DO jj=nldj, nlej
167               DO ji=nldi, nlei
168                  pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj)
169               END DO
170            END DO
171            CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. )
172         END DO
173         
174         ! temperature:
175         pts(:,:,:,jp_tem) = 10._wp
176         ! salinity: 
177         DO jk=1, jpkm1
178            pts(:,:,jk,jp_sal) = gphit(:,:)
179         END DO
180         ! velocities:
181         pv(:,:,:) = 0.
182         !           
183      CASE(4)    ! geostrophic zonal pulse
184   
185         DO jj=1, jpj
186            DO ji=1, jpi
187               IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN
188                  zdu = rn_uzonal
189               ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN
190                  zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. )
191               ELSE
192                  zdu = 0.
193               END IF
194               IF ( ABS(gphit(ji,jj)) <= zjety ) THEN
195                  pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav
196                  pu(ji,jj,:) = zdu
197                  pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1.
198               ELSE
199                  pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav 
200                  pu(ji,jj,:) = 0.
201                  pts(ji,jj,:,jp_sal) = 1.
202               END IF
203            END DO
204         END DO
205         
206         ! temperature:
207         pts(:,:,:,jp_tem) = 10._wp * ptmask(:,:,:)       
208         pv(:,:,:) = 0.
209         
210         
211       CASE(5)    ! vortex
212                  !
213         zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 )
214         zumax = rn_vtxmax * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic
215         zlambda = SQRT(2._wp)*rn_lambda       ! Horizontal scale in meters
216         zn2 = 3.e-3**2
217         zH = 0.5_wp * 5000._wp
218         !
219         zr_lambda2 = 1._wp / zlambda**2
220         zP0 = rho0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
221         !
222         DO jj=1, jpj
223            DO ji=1, jpi
224               zx = glamt(ji,jj) * 1.e3
225               zy = gphit(ji,jj) * 1.e3
226               ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y)
227               zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy
228               ! Sea level:
229               pssh(ji,jj) = 0.
230               DO jl=1,5
231                  zdt = pssh(ji,jj)
232                  zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z)
233                  zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)
234                  pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g)
235               END DO
236               ! temperature:
237               DO jk=1,jpk
238                  zdt =  pdept(ji,jj,jk) 
239                  zrho1 = rho0 * (1._wp + zn2*zdt/grav)
240                  IF (zdt < zH) THEN
241                     zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH))   ! F'(z)
242                     zrho1 = zrho1 - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)
243                  ENDIF
244                  !               pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
245                  pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
246               END DO
247            END DO
248         END DO
249         !
250         ! salinity: 
251         pts(:,:,:,jp_sal) = 35._wp * ptmask(:,:,:) 
252         !
253         ! velocities:
254         za = 2._wp * zP0 / zlambda**2
255         DO jj=1, jpj
256            DO ji=1, jpim1
257               zx = glamu(ji,jj) * 1.e3
258               zy = gphiu(ji,jj) * 1.e3
259               DO jk=1, jpk
260                  zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk))
261                  IF (zdu < zH) THEN
262                     zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH))
263                     zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal
264                     pu(ji,jj,jk) = - zf / ( rho0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk)
265                  ELSE
266                     pu(ji,jj,jk) = 0._wp
267                  ENDIF
268               END DO
269            END DO
270         END DO
271         !
272         DO jj=1, jpjm1
273            DO ji=1, jpi
274               zx = glamv(ji,jj) * 1.e3
275               zy = gphiv(ji,jj) * 1.e3
276               DO jk=1, jpk
277                  zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk))
278                  IF (zdv < zH) THEN
279                     zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH))
280                     zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2)
281                     pv(ji,jj,jk) = zf / ( rho0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk)
282                  ELSE
283                     pv(ji,jj,jk) = 0._wp
284                  ENDIF
285               END DO
286            END DO
287         END DO
288         !           
289      END SELECT
290
291      IF (ln_sshnoise) THEN
292         CALL RANDOM_NUMBER(zrandom)
293         pssh(:,:) = pssh(:,:) + ( 0.1  * zrandom(:,:) - 0.05 )
294      END IF
295      CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. )
296      CALL lbc_lnk(  'usrdef_istate', pts, 'T',  1. )
297      CALL lbc_lnk(   'usrdef_istate', pu, 'U', -1. )
298      CALL lbc_lnk(   'usrdef_istate', pv, 'V', -1. )
299
300   END SUBROUTINE usr_def_istate
301
302   !!======================================================================
303END MODULE usrdef_istate
Note: See TracBrowser for help on using the repository browser.