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/2020/dev_r12558_HPC-08_epico_Extra_Halo/tests/CANAL/MY_SRC – NEMO

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/tests/CANAL/MY_SRC/usrdef_istate.F90 @ 12939

Last change on this file since 12939 was 12939, checked in by smasson, 4 years ago

Extra_Halo: update with trunk@12933, see #2366
  • Property svn:keywords set to Id
File size: 11.8 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   !! * Substitutions
31#  include "do_loop_substitute.h90"
32   !!----------------------------------------------------------------------
33   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
34   !! $Id$
35   !! Software governed by the CeCILL license (see ./LICENSE)
36   !!----------------------------------------------------------------------
37CONTAINS
38 
39   SUBROUTINE usr_def_istate( pdept, ptmask, pts, pu, pv, pssh )
40      !!----------------------------------------------------------------------
41      !!                   ***  ROUTINE usr_def_istate  ***
42      !!
43      !! ** Purpose :   Initialization of the dynamics and tracers
44      !!                Here CANAL configuration
45      !!
46      !! ** Method  :   Set a gaussian anomaly of pressure and associated
47      !!                geostrophic velocities
48      !!----------------------------------------------------------------------
49      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pdept   ! depth of t-point               [m]
50      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   ptmask  ! t-point ocean mask             [m]
51      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pts     ! T & S fields      [Celsius ; g/kg]
52      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pu      ! i-component of the velocity  [m/s]
53      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pv      ! j-component of the velocity  [m/s]
54      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(  out) ::   pssh    ! sea-surface height
55      !
56      INTEGER  :: ji, jj, jk, jl  ! dummy loop indices
57      REAL(wp) :: zx, zy, zP0, zumax, zlambda, zr_lambda2, zn2, zf0, zH, zrho1, za, zf, zdzF
58      REAL(wp) :: zpsurf, zdyPs, zdxPs
59      REAL(wp) :: zdt, zdu, zdv
60      REAL(wp) :: zjetx, zjety, zbeta
61      REAL(wp), DIMENSION(jpi,jpj)  ::   zrandom
62      !!----------------------------------------------------------------------
63      !
64      IF(lwp) WRITE(numout,*)
65      IF(lwp) WRITE(numout,*) 'usr_def_istate : CANAL configuration, analytical definition of initial state'
66      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   '
67      !
68      IF (ln_sshnoise) CALL RANDOM_NUMBER(zrandom)
69      zjetx = ABS(rn_ujetszx)/2.
70      zjety = ABS(rn_ujetszy)/2.
71      !
72      SELECT CASE(nn_initcase)
73      CASE(0)    ! rest
74         
75         ! sea level:
76         pssh(:,:) = 0.
77         ! temperature:
78         pts(:,:,:,jp_tem) = 10._wp
79         ! salinity: 
80         pts(:,:,:,jp_sal) = 35._wp
81         ! velocities:
82         pu(:,:,:) = 0.
83         pv(:,:,:) = 0.
84         
85      CASE(1)    ! geostrophic zonal jet from -zjety to +zjety
86
87         ! sea level:
88         SELECT CASE( nn_fcase )
89         CASE(0)    ! f = f0
90            ! sea level: ssh = - fuy / g
91            WHERE( ABS(gphit) <= zjety )
92               pssh(:,:) = - ff_t(:,:) * rn_uzonal * gphit(:,:) * 1.e3 / grav
93            ELSEWHERE
94               pssh(:,:) = - ff_t(:,:) * rn_uzonal * SIGN(zjety, gphit(:,:)) * 1.e3 / grav
95            END WHERE
96         CASE(1)    ! f = f0 + beta*y
97            ! sea level: ssh = - u / g * ( fy + 0.5 * beta * y^2 )
98            zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra
99            WHERE( ABS(gphit) <= zjety )
100               pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
101            ELSEWHERE
102               pssh(:,:) = - rn_uzonal / grav * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3   &
103                  &                             + 0.5 * zbeta * zjety * zjety * 1.e6 )
104            END WHERE
105         END SELECT
106         ! temperature:
107         pts(:,:,:,jp_tem) = 10._wp
108         ! salinity: 
109         pts(:,:,jpk,jp_sal) = 0.
110         DO jk=1, jpkm1
111            pts(:,:,jk,jp_sal) = gphit(:,:)
112         END DO
113         ! velocities:
114         pu(:,:,:) = 0.
115         DO jk=1, jpkm1
116            WHERE( ABS(gphit) <= zjety ) pu(:,:,jk) = rn_uzonal
117         END DO
118         pv(:,:,:) = 0.
119         !                 
120      CASE(2)    ! geostrophic zonal current shear
121     
122         ! sea level:
123         SELECT CASE( nn_fcase )
124         CASE(0)    ! f = f0
125            ! sea level: ssh = - fuy / g
126            WHERE( ABS(gphit) <= zjety )
127               pssh(:,:) = - ff_t(:,:) * rn_uzonal * ABS(gphit(:,:)) * 1.e3 / grav
128            ELSEWHERE
129               pssh(:,:) = - ff_t(:,:) * rn_uzonal * zjety * 1.e3 / grav
130            END WHERE
131         CASE(1)    ! f = f0 + beta*y
132            ! sea level: ssh = - u / g * ( fy + 0.5 * beta * y^2 )
133            zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra
134            WHERE( ABS(gphit) <= zjety )
135               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   &
136                  &        * ( ff_t(:,:) * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
137            ELSEWHERE
138               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   &
139                  &        * ( ff_t(:,:) * SIGN(zjety, gphit(:,:)) * 1.e3 + 0.5 * zbeta * zjety * zjety * 1.e6 )
140            END WHERE
141         END SELECT
142         ! temperature:
143         pts(:,:,:,jp_tem) = 10._wp
144         ! salinity: 
145         pts(:,:,:,jp_sal) = 2.
146         DO jk=1, jpkm1
147            WHERE( ABS(gphiv) <= zjety ) pts(:,:,jk,jp_sal) = 2. + SIGN(1.,gphiv(:,:))
148         END DO
149         ! velocities:
150         pu(:,:,:) = 0.
151         DO jk=1, jpkm1
152            WHERE( ABS(gphiv) <= zjety ) pu(:,:,jk) = SIGN(rn_uzonal,gphit(:,:))*SIGN(1.,rn_uzonal)
153            WHERE( ABS(gphiv) == 0.    ) pu(:,:,jk) = 0. 
154         END DO
155         pv(:,:,:) = 0.
156         !                 
157      CASE(3)    ! gaussian zonal currant
158
159         ! zonal current
160         DO jk=1, jpkm1
161            ! gphit and lambda are both in km
162            pu(:,:,jk) = rn_uzonal * EXP( - 0.5 * gphit(:,:)**2 / rn_lambda**2 )
163         END DO
164         
165         ! sea level:
166         pssh(:,1) = - ff_t(:,1) / grav * pu(:,1,1) * e2t(:,1)
167         DO jl=1, jpnj
168            DO_2D_00_00
169               pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj)
170            END_2D
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_2D_11_11
186            IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN
187               zdu = rn_uzonal
188            ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN
189               zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. )
190            ELSE
191               zdu = 0.
192            END IF
193            IF ( ABS(gphit(ji,jj)) <= zjety ) THEN
194               pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav
195               pu(ji,jj,:) = zdu
196               pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1.
197            ELSE
198               pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav 
199               pu(ji,jj,:) = 0.
200               pts(ji,jj,:,jp_sal) = 1.
201            END IF
202         END_2D
203         
204         ! temperature:
205         pts(:,:,:,jp_tem) = 10._wp * ptmask(:,:,:)       
206         pv(:,:,:) = 0.
207         
208       CASE(5)    ! vortex
209                  !
210         zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 )
211         zumax = rn_vtxmax * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic
212         zlambda = SQRT(2._wp)*rn_lambda       ! Horizontal scale in meters
213         zn2 = 3.e-3**2
214         zH = 0.5_wp * 5000._wp
215         !
216         zr_lambda2 = 1._wp / zlambda**2
217         zP0 = rho0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
218         !
219         DO_2D_11_11
220            zx = glamt(ji,jj) * 1.e3
221            zy = gphit(ji,jj) * 1.e3
222            ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y)
223            zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy
224            ! Sea level:
225            pssh(ji,jj) = 0.
226            DO jl=1,5
227               zdt = pssh(ji,jj)
228               zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z)
229               zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)
230               pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g)
231            END DO
232            ! temperature:
233            DO jk=1,jpk
234               zdt =  pdept(ji,jj,jk) 
235               zrho1 = rho0 * (1._wp + zn2*zdt/grav)
236               IF (zdt < zH) THEN
237                  zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH))   ! F'(z)
238                  zrho1 = zrho1 - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)
239               ENDIF
240               !               pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
241               pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
242            END DO
243         END_2D
244         !
245         ! salinity: 
246         pts(:,:,:,jp_sal) = 35._wp * ptmask(:,:,:) 
247         !
248         ! velocities:
249         za = 2._wp * zP0 / zlambda**2
250         DO_2D_00_00
251            zx = glamu(ji,jj) * 1.e3
252            zy = gphiu(ji,jj) * 1.e3
253            DO jk=1, jpk
254               zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk))
255               IF (zdu < zH) THEN
256                  zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH))
257                  zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal
258                  pu(ji,jj,jk) = - zf / ( rho0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk)
259               ELSE
260                  pu(ji,jj,jk) = 0._wp
261               ENDIF
262            END DO
263         END_2D
264         !
265         DO_2D_00_00
266            zx = glamv(ji,jj) * 1.e3
267            zy = gphiv(ji,jj) * 1.e3
268            DO jk=1, jpk
269               zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk))
270               IF (zdv < zH) THEN
271                  zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH))
272                  zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2)
273                  pv(ji,jj,jk) = zf / ( rho0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk)
274               ELSE
275                  pv(ji,jj,jk) = 0._wp
276               ENDIF
277            END DO
278         END_2D
279         !           
280      END SELECT
281     
282      IF (ln_sshnoise) THEN
283         CALL RANDOM_NUMBER(zrandom)
284         pssh(:,:) = pssh(:,:) + ( 0.1  * zrandom(:,:) - 0.05 )
285      END IF
286      CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. )
287      CALL lbc_lnk( 'usrdef_istate', pts , 'T',  1. )
288      CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. )
289
290   END SUBROUTINE usr_def_istate
291
292   !!======================================================================
293END MODULE usrdef_istate
Note: See TracBrowser for help on using the repository browser.