source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/tests/CANAL/MY_SRC/usrdef_istate.F90 @ 12751

Last change on this file since 12751 was 12751, checked in by clem, 9 months ago

update CANAL for students experiments

  • Property svn:keywords set to Id
File size: 12.6 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 jj=nldj, nlej
187               DO ji=nldi, nlei
188                  pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj)
189               END DO
190            END DO
191            CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. )
192         END DO
193         
194         ! temperature:
195         pts(:,:,:,jp_tem) = 10._wp
196         ! salinity: 
197         DO jk=1, jpkm1
198            pts(:,:,jk,jp_sal) = pssh(:,:)
199         END DO
200         ! velocities:
201         pv(:,:,:) = 0.
202         !           
203      CASE(4)    ! geostrophic zonal pulse
204   
205         DO jj=1, jpj
206            DO ji=1, jpi
207               IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN
208                  zdu = rn_uzonal
209               ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN
210                  zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. )
211               ELSE
212                  zdu = 0.
213               END IF
214               IF ( ABS(gphit(ji,jj)) <= zjety ) THEN
215                  pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav
216                  pu(ji,jj,:) = zdu
217                  pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1.
218               ELSE
219                  pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav 
220                  pu(ji,jj,:) = 0.
221                  pts(ji,jj,:,jp_sal) = 1.
222               END IF
223            END DO
224         END DO
225         
226         ! temperature:
227         pts(:,:,:,jp_tem) = 10._wp * ptmask(:,:,:)       
228         pv(:,:,:) = 0.
229         
230         
231       CASE(5)    ! vortex
232                  !
233         zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 )
234         zumax = rn_vtxmax * SIGN(1._wp, zf0) ! Here Anticyclonic: set zumax=-1 for cyclonic
235         zlambda = SQRT(2._wp)*rn_lambda*1.e3       ! Horizontal scale in meters
236         zn2 = 3.e-3**2
237         zH = 0.5_wp * 5000._wp
238         !
239         zr_lambda2 = 1._wp / zlambda**2
240         zP0 = rau0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
241         !
242         DO jj=1, jpj
243            DO ji=1, jpi
244               zx = glamt(ji,jj) * 1.e3
245               zy = gphit(ji,jj) * 1.e3
246               ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y)
247               zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rau0 * ff_t(ji,jj) * rn_uzonal * zy
248               ! Sea level:
249               pssh(ji,jj) = 0.
250               DO jl=1,5
251                  zdt = pssh(ji,jj)
252                  zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z)
253                  zrho1 = rau0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)
254                  pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g)
255               END DO
256               ! temperature:
257               DO jk=1,jpk
258                  zdt =  pdept(ji,jj,jk) 
259                  zrho1 = rau0 * (1._wp + zn2*zdt/grav)
260                  IF (zdt < zH) THEN
261                     zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH))   ! F'(z)
262                     zrho1 = zrho1 - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y)
263                  ENDIF
264                  !               pts(ji,jj,jk,jp_tem) = (20._wp + (rau0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
265                  pts(ji,jj,jk,jp_tem) = (10._wp + (rau0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk)
266               END DO
267            END DO
268         END DO
269         !
270         ! salinity: 
271         pts(:,:,:,jp_sal) = 35._wp * ptmask(:,:,:) 
272         !
273         ! velocities:
274         za = 2._wp * zP0 / zlambda**2
275         DO jj = 2, jpjm1
276            DO ji = 2, jpim1
277               zx = glamu(ji,jj) * 1.e3
278               zy = gphiu(ji,jj) * 1.e3
279               DO jk=1, jpk
280                  zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk))
281                  IF (zdu < zH) THEN
282                     zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH))
283                     zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rau0 * ff_t(ji,jj) * rn_uzonal
284                     pu(ji,jj,jk) = - zf / ( rau0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk)
285                  ELSE
286                     pu(ji,jj,jk) = 0._wp
287                  ENDIF
288               END DO
289            END DO
290         END DO
291         !
292         DO jj = 2, jpjm1
293            DO ji = 2, jpim1
294               zx = glamv(ji,jj) * 1.e3
295               zy = gphiv(ji,jj) * 1.e3
296               DO jk=1, jpk
297                  zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk))
298                  IF (zdv < zH) THEN
299                     zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH))
300                     zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2)
301                     pv(ji,jj,jk) = zf / ( rau0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk)
302                  ELSE
303                     pv(ji,jj,jk) = 0._wp
304                  ENDIF
305               END DO
306            END DO
307         END DO
308         !           
309         CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. )
310
311      END SELECT
312
313      IF (ln_sshnoise) THEN
314         CALL RANDOM_SEED()
315         CALL RANDOM_NUMBER(zrandom)
316         pssh(:,:) = pssh(:,:) + ( 0.1  * zrandom(:,:) - 0.05 )
317      END IF
318 
319   END SUBROUTINE usr_def_istate
320
321   !!======================================================================
322END MODULE usrdef_istate
Note: See TracBrowser for help on using the repository browser.