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/temporary_r4_trunk/tests/CANAL/MY_SRC – NEMO

source: NEMO/branches/2020/temporary_r4_trunk/tests/CANAL/MY_SRC/usrdef_istate.F90 @ 13469

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

r4_trunk: first change of DO loops for routines to be merged, see #2523

  • 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      zjetx = ABS(rn_ujetszx)/2.
67      zjety = ABS(rn_ujetszy)/2.
68      !
69      zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 )
70      !
71      SELECT CASE(nn_initcase)
72
73      CASE(-1)    ! stratif at rest
74
75         ! sea level:
76         pssh(:,:) = 0.
77         ! temperature:
78         pts(:,:,1,jp_tem) = 25. !!30._wp
79         pts(:,:,2:jpk,jp_tem) = 22. !!24._wp
80         ! salinity: 
81         pts(:,:,:,jp_sal) = 35._wp
82         ! velocities:
83         pu(:,:,:) = 0.
84         pv(:,:,:) = 0.
85
86      CASE(0)    ! rest
87         
88         ! sea level:
89         pssh(:,:) = 0.
90         ! temperature:
91         pts(:,:,:,jp_tem) = 10._wp
92         ! salinity: 
93         pts(:,:,:,jp_sal) = 35._wp
94         ! velocities:
95         pu(:,:,:) = 0.
96         pv(:,:,:) = 0.
97         
98      CASE(1)    ! geostrophic zonal jet from -zjety to +zjety
99
100         ! sea level:
101         SELECT CASE( nn_fcase )
102         CASE(0)    ! f = f0
103            ! sea level: ssh = - fuy / g
104            WHERE( ABS(gphit) <= zjety )
105               pssh(:,:) = - ff_t(:,:) * rn_uzonal * gphit(:,:) * 1.e3 / grav
106            ELSEWHERE
107               pssh(:,:) = - ff_t(:,:) * rn_uzonal * SIGN(zjety, gphit(:,:)) * 1.e3 / grav
108            END WHERE
109         CASE(1)    ! f = f0 + beta*y
110            ! sea level: ssh = - u / g * ( fy + 0.5 * beta * y^2 )
111            zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra
112            WHERE( ABS(gphit) <= zjety )
113               pssh(:,:) = - rn_uzonal / grav * ( zf0 * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
114            ELSEWHERE
115               pssh(:,:) = - rn_uzonal / grav * ( zf0 * SIGN(zjety, gphit(:,:)) * 1.e3   &
116                  &                             + 0.5 * zbeta * zjety * zjety * 1.e6 )
117            END WHERE
118         END SELECT
119         ! temperature:
120         pts(:,:,:,jp_tem) = 10._wp
121         ! salinity: 
122         pts(:,:,jpk,jp_sal) = 0.
123         DO jk=1, jpkm1
124            WHERE( ABS(gphit) <= zjety )
125!!$            WHERE( ABS(gphit) <= zjety*0.5 .AND. ABS(glamt) <= zjety*0.5 ) ! for a square of salt
126               pts(:,:,jk,jp_sal) = 35.
127            ELSEWHERE
128               pts(:,:,jk,jp_sal) = 30.
129            END WHERE                   
130         END DO
131         ! velocities:
132         pu(:,:,:) = 0.
133         DO jk=1, jpkm1
134            WHERE( ABS(gphit) <= zjety ) pu(:,:,jk) = rn_uzonal
135         END DO
136         pv(:,:,:) = 0.
137         !                 
138      CASE(2)    ! geostrophic zonal current shear
139     
140         ! sea level:
141         SELECT CASE( nn_fcase )
142         CASE(0)    ! f = f0
143            ! sea level: ssh = - fuy / g
144            WHERE( ABS(gphit) <= zjety )
145               pssh(:,:) = - ff_t(:,:) * rn_uzonal * ABS(gphit(:,:)) * 1.e3 / grav
146            ELSEWHERE
147               pssh(:,:) = - ff_t(:,:) * rn_uzonal * zjety * 1.e3 / grav
148            END WHERE
149         CASE(1)    ! f = f0 + beta*y
150            ! sea level: ssh = - u / g * ( fy + 0.5 * beta * y^2 )
151            zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra
152            WHERE( ABS(gphit) <= zjety )
153               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   &
154                  &        * ( zf0 * gphit(:,:) * 1.e3 + 0.5 * zbeta * gphit(:,:) * gphit(:,:) * 1.e6 )
155            ELSEWHERE
156               pssh(:,:) = - SIGN(rn_uzonal, gphit(:,:)) / grav   &
157                  &        * ( zf0 * SIGN(zjety, gphit(:,:)) * 1.e3 + 0.5 * zbeta * zjety * zjety * 1.e6 )
158            END WHERE
159         END SELECT
160         ! temperature:
161         pts(:,:,:,jp_tem) = 10._wp
162         ! salinity: 
163         pts(:,:,:,jp_sal) = 30.
164         DO jk=1, jpkm1
165            WHERE( ABS(gphiv) <= zjety ) pts(:,:,jk,jp_sal) = 30. + SIGN(1.,gphiv(:,:))
166         END DO
167         ! velocities:
168         pu(:,:,:) = 0.
169         DO jk=1, jpkm1
170            WHERE( ABS(gphiv) <= zjety ) pu(:,:,jk) = SIGN(rn_uzonal,gphit(:,:))*SIGN(1.,rn_uzonal)
171            WHERE( ABS(gphiv) == 0.    ) pu(:,:,jk) = 0. 
172         END DO
173         pv(:,:,:) = 0.
174         !                 
175      CASE(3)    ! gaussian zonal currant
176
177         ! zonal current
178         DO jk=1, jpkm1
179            ! gphit and lambda are both in km
180            pu(:,:,jk) = rn_uzonal * EXP( - 0.5 * gphit(:,:)**2 / rn_lambda**2 )
181         END DO
182         
183         ! sea level:
184         pssh(:,1) = - ff_t(:,1) / grav * pu(:,1,1) * e2t(:,1)
185         DO jl=1, jpnj
186            DO_2D_nldj1_nldi1
187               pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj)
188            END_2D
189            CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. )
190         END DO
191         
192         ! temperature:
193         pts(:,:,:,jp_tem) = 10._wp
194         ! salinity: 
195         DO jk=1, jpkm1
196            pts(:,:,jk,jp_sal) = pssh(:,:)
197         END DO
198         ! velocities:
199         pv(:,:,:) = 0.
200         !           
201      CASE(4)    ! geostrophic zonal pulse
202   
203         DO_2D_11_11
204            IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN
205               zdu = rn_uzonal
206            ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN
207               zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. )
208            ELSE
209               zdu = 0.
210            END IF
211            IF ( ABS(gphit(ji,jj)) <= zjety ) THEN
212               pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav
213               pu(ji,jj,:) = zdu
214               pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1.
215            ELSE
216               pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav 
217               pu(ji,jj,:) = 0.
218               pts(ji,jj,:,jp_sal) = 1.
219            END IF
220         END_2D
221         
222         ! temperature:
223         pts(:,:,:,jp_tem) = 10._wp * ptmask(:,:,:)       
224         pv(:,:,:) = 0.
225         
226         
227       CASE(5)    ! vortex
228                  !
229         zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 )
230         zumax = rn_vtxmax * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic
231         zlambda = SQRT(2._wp)*rn_lambda*1.e3       ! Horizontal scale in meters
232         zn2 = 3.e-3**2
233         zH = 0.5_wp * 5000._wp
234         !
235         zr_lambda2 = 1._wp / zlambda**2
236         zP0 = rau0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
237         !
238         DO_2D_11_11
239            zx = glamt(ji,jj) * 1.e3
240            zy = gphit(ji,jj) * 1.e3
241            ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y)
242            zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rau0 * ff_t(ji,jj) * rn_uzonal * zy
243            ! Sea level:
244            pssh(ji,jj) = 0.
245            DO jl=1,5
246               zdt = pssh(ji,jj)
247               zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z)
248               zrho1 = rau0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)
249               pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g)
250            END DO
251            ! temperature:
252            DO jk=1,jpk
253               zdt =  pdept(ji,jj,jk) 
254               zrho1 = rau0 * (1._wp + zn2*zdt/grav)
255               IF (zdt < zH) THEN
256                  zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH))   ! F'(z)
257                  zrho1 = zrho1 - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)
258               ENDIF
259               !               pts(ji,jj,jk,jp_tem) = (20._wp + (rau0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
260               pts(ji,jj,jk,jp_tem) = (10._wp + (rau0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
261            END DO
262         END_2D
263         !
264         ! salinity: 
265         pts(:,:,:,jp_sal) = 35._wp * ptmask(:,:,:) 
266         !
267         ! velocities:
268         za = 2._wp * zP0 / zlambda**2
269         DO_2D_00_00
270            zx = glamu(ji,jj) * 1.e3
271            zy = gphiu(ji,jj) * 1.e3
272            DO jk=1, jpk
273               zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk))
274               IF (zdu < zH) THEN
275                  zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH))
276                  zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rau0 * ff_t(ji,jj) * rn_uzonal
277                  pu(ji,jj,jk) = - zf / ( rau0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk)
278               ELSE
279                  pu(ji,jj,jk) = 0._wp
280               ENDIF
281            END DO
282         END_2D
283         !
284         DO_2D_00_00
285            zx = glamv(ji,jj) * 1.e3
286            zy = gphiv(ji,jj) * 1.e3
287            DO jk=1, jpk
288               zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk))
289               IF (zdv < zH) THEN
290                  zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH))
291                  zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2)
292                  pv(ji,jj,jk) = zf / ( rau0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk)
293               ELSE
294                  pv(ji,jj,jk) = 0._wp
295               ENDIF
296            END DO
297         END_2D
298         !           
299         CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. )
300
301      END SELECT
302
303      IF (ln_sshnoise) THEN
304         CALL RANDOM_SEED()
305         CALL RANDOM_NUMBER(zrandom)
306         pssh(:,:) = pssh(:,:) + ( 0.1  * zrandom(:,:) - 0.05 )
307      END IF
308 
309   END SUBROUTINE usr_def_istate
310
311   !!======================================================================
312END MODULE usrdef_istate
Note: See TracBrowser for help on using the repository browser.