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.
domvvl.F90 in trunk/NEMO/OPA_SRC/DOM – NEMO

source: trunk/NEMO/OPA_SRC/DOM/domvvl.F90 @ 627

Last change on this file since 627 was 592, checked in by opalod, 17 years ago

nemo_v2_update_001 : CT : - add non linear free surface (variable volume) with new cpp key key_vvl
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.9 KB
Line 
1MODULE domvvl
2   !!======================================================================
3   !!                       ***  MODULE domvvl   ***
4   !! Ocean :
5   !!======================================================================
6   !! History :   9.0  !  06-06  (B. Levier, L. Marie)  original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   'key_vvl'                              variable volume
11   !!----------------------------------------------------------------------
12   !!----------------------------------------------------------------------
13   !!   dom_vvl         : defined scale factors & depths at each time step
14   !!   dom_vvl_ini     : defined coefficients to distribute ssh on each layers
15   !!   dom_vvl_ssh     : defined the ocean sea level at each time step
16   !!----------------------------------------------------------------------
17   !! * Modules used
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE in_out_manager  ! I/O manager
21   USE lib_mpp         ! distributed memory computing library
22   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
23   USE dynspg_oce      ! surface pressure gradient variables
24   USE ocesbc          ! ocean surface boundary condition
25   USE phycst          ! physical constants
26
27   IMPLICIT NONE
28   PRIVATE
29
30   !! * Routine accessibility
31   PUBLIC dom_vvl_ini    ! called by dom_init.F90
32   PUBLIC dom_vvl_ssh    ! called by trazdf.F90
33   PUBLIC dom_vvl        ! called by istate.F90 and step.F90
34
35   !! * Module variables
36   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !:
37      mut, muu, muv, muf                            !:
38
39   REAL(wp), DIMENSION(jpk) ::   r2dt               ! vertical profile time-step, = 2 rdttra
40      !                                             ! except at nit000 (=rdttra) if neuler=0
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !!   OPA 9.0 , LOCEAN-IPSL (2005)
47   !! $Header$
48   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50
51CONTAINS       
52
53#if defined key_vvl
54
55   SUBROUTINE dom_vvl_ini
56      !!----------------------------------------------------------------------
57      !!                ***  ROUTINE dom_vvl_ini  ***
58      !!                   
59      !! ** Purpose :  compute coefficients muX at T-U-V-F points to spread
60      !!               ssh over the whole water column (scale factors)
61      !!
62      !!----------------------------------------------------------------------
63      INTEGER  ::   ji, jj, jk, zpk
64      REAL(wp), DIMENSION(jpi,jpj) ::   zmut, zmuu, zmuv, zmuf   ! 2D workspace
65      !!----------------------------------------------------------------------
66
67      IF(lwp)   THEN
68         WRITE(numout,*)
69         WRITE(numout,*) 'dom_vvl_ini : Variable volume activated'
70         WRITE(numout,*) '~~~~~~~~~~~   compute coef. used to spread ssh over each layers'
71      ENDIF
72
73      IF( ln_zps )  CALL ctl_stop( 'dom_vvl_ini : option ln_zps is incompatible with variable volume option key_vvl')
74
75#if defined key_zco  ||  defined key_dynspg_rl
76      CALL ctl_stop( 'dom_vvl_ini : options key_zco/key_dynspg_rl are incompatible with variable volume option key_vvl')
77#endif
78
79      fsvdept (:,:,:) = gdept (:,:,:)
80      fsvdepw (:,:,:) = gdepw (:,:,:)
81      fsvde3w (:,:,:) = gdep3w(:,:,:)
82      fsve3t  (:,:,:) = e3t   (:,:,:)
83      fsve3u  (:,:,:) = e3u   (:,:,:)
84      fsve3v  (:,:,:) = e3v   (:,:,:)
85      fsve3f  (:,:,:) = e3f   (:,:,:)
86      fsve3w  (:,:,:) = e3w   (:,:,:)
87      fsve3uw (:,:,:) = e3uw  (:,:,:)
88      fsve3vw (:,:,:) = e3vw  (:,:,:)
89
90      ! mu computation
91      zmut(:,:)   = 0.e0
92      zmuu(:,:)   = 0.e0
93      zmuv(:,:)   = 0.e0
94      zmuf(:,:)   = 0.e0
95      mut (:,:,:) = 0.e0
96      muu (:,:,:) = 0.e0
97      muv (:,:,:) = 0.e0
98      muf (:,:,:) = 0.e0
99
100      DO jj = 1, jpj
101         DO ji = 1, jpi
102            zpk = mbathy(ji,jj) - 1
103            DO jk = 1, zpk
104               zmut(ji,jj) = zmut(ji,jj) + fsve3t(ji,jj,jk) * SUM( fsve3t(ji,jj,jk:zpk) )
105               zmuu(ji,jj) = zmuu(ji,jj) + fsve3u(ji,jj,jk) * SUM( fsve3u(ji,jj,jk:zpk) )
106               zmuv(ji,jj) = zmuv(ji,jj) + fsve3v(ji,jj,jk) * SUM( fsve3v(ji,jj,jk:zpk) )
107               zmuf(ji,jj) = zmuf(ji,jj) + fsve3f(ji,jj,jk) * SUM( fsve3f(ji,jj,jk:zpk) )
108            END DO
109         END DO
110      END DO
111
112      DO jj = 1, jpj
113         DO ji = 1, jpi
114            zpk = mbathy(ji,jj) - 1
115            DO jk = 1, zpk
116               mut(ji,jj,jk) = SUM( fsve3t(ji,jj,jk:zpk) ) / zmut(ji,jj)
117               muu(ji,jj,jk) = SUM( fsve3u(ji,jj,jk:zpk) ) / zmuu(ji,jj)
118               muv(ji,jj,jk) = SUM( fsve3v(ji,jj,jk:zpk) ) / zmuv(ji,jj)
119               muf(ji,jj,jk) = SUM( fsve3f(ji,jj,jk:zpk) ) / zmuf(ji,jj)
120            END DO
121         END DO
122      END DO
123
124
125   END SUBROUTINE dom_vvl_ini
126
127
128
129   SUBROUTINE dom_vvl
130      !!----------------------------------------------------------------------
131      !!                ***  ROUTINE dom_vvl  ***
132      !!                   
133      !! ** Purpose :  compute ssh at U-V-F points, T-W scale factors and local
134      !!               depths at each time step.
135      !!----------------------------------------------------------------------
136      !! * Local declarations
137      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
138      REAL(wp), DIMENSION(jpi,jpj) :: zsshf    ! 2D workspace
139      !!----------------------------------------------------------------------
140
141      ! Sea Surface Height at u-v-fpoints
142      DO jj = 1, jpjm1
143         DO ji = 1,jpim1
144            sshu(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )  &
145               &          * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
146               &            + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
147            !
148            sshv(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )  &
149               &          * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
150               &            + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
151            !
152            zsshf(ji,jj) = 0.25 * fmask(ji,jj,1)                 &
153               &           * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   &
154               &             + sshn(ji+1,jj) + sshn(ji+1,jj+1) )
155         END DO
156      END DO
157
158      ! Boundaries conditions
159      CALL lbc_lnk( sshu,  'U', 1. )
160      CALL lbc_lnk( sshv,  'V', 1. )
161      CALL lbc_lnk( zsshf, 'F', 1. )
162
163      ! Scale factors at T levels
164      DO jk = 1, jpkm1
165         fse3t(:,:,jk) = fsve3t(:,:,jk) * ( 1 + sshn(:,:)  * mut(:,:,jk) )
166         fse3u(:,:,jk) = fsve3u(:,:,jk) * ( 1 + sshu(:,:)  * muu(:,:,jk) )
167         fse3v(:,:,jk) = fsve3v(:,:,jk) * ( 1 + sshv(:,:)  * muv(:,:,jk) )
168         fse3f(:,:,jk) = fsve3f(:,:,jk) * ( 1 + zsshf(:,:) * muf(:,:,jk) )
169      END DO
170
171      ! Scale factors at W levels
172      fse3w (:,:,1) = fse3t(:,:,1)
173      fse3uw(:,:,1) = fse3u(:,:,1)
174      fse3vw(:,:,1) = fse3v(:,:,1)
175      DO jk = 2, jpk
176         fse3w (:,:,jk) = 0.5 * ( fse3t(:,:,jk-1) + fse3t(:,:,jk) )
177         fse3uw(:,:,jk) = 0.5 * ( fse3u(:,:,jk-1) + fse3u(:,:,jk) )
178         fse3vw(:,:,jk) = 0.5 * ( fse3v(:,:,jk-1) + fse3v(:,:,jk) )
179      END DO
180
181      ! T and W points depth
182      fsdept(:,:,1) = 0.5 * fse3w(:,:,1)
183      fsdepw(:,:,1) = 0.e0
184      fsde3w(:,:,1) = fsdept(:,:,1) - sshn(:,:)
185      DO jk = 2, jpk
186         fsdept(:,:,jk) = fsdept(:,:,jk-1) + fse3w(:,:,jk)
187         fsdepw(:,:,jk) = fsdepw(:,:,jk-1) + fse3t(:,:,jk-1)
188         fsde3w(:,:,jk) = fsdept(:,:,jk  ) - sshn(:,:)
189      END DO
190
191      ! Local depth or Inverse of the local depth of the water column at u- and v-points
192      ! ------------------------------
193
194      ! Ocean depth at U- and V-points
195      hu(:,:) = 0.e0    ;    hv(:,:) = 0.e0
196
197      DO jk = 1, jpk
198         hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk)
199         hv(:,:) = hv(:,:) + fse3v(:,:,jk) * vmask(:,:,jk)
200      END DO
201
202
203      ! Inverse of the local depth
204      hur(:,:) = fse3u(:,:,1)             ! Lower bound : thickness of the first model level
205      hvr(:,:) = fse3v(:,:,1)
206     
207      DO jk = 2, jpk                      ! Sum of the vertical scale factors
208         hur(:,:) = hur(:,:) + fse3u(:,:,jk) * umask(:,:,jk)
209         hvr(:,:) = hvr(:,:) + fse3v(:,:,jk) * vmask(:,:,jk)
210      END DO
211
212      ! Compute and mask the inverse of the local depth
213      hur(:,:) = 1. / hur(:,:) * umask(:,:,1)
214      hvr(:,:) = 1. / hvr(:,:) * vmask(:,:,1)
215
216   END SUBROUTINE dom_vvl
217
218
219
220   SUBROUTINE dom_vvl_ssh( kt ) 
221      !!----------------------------------------------------------------------
222      !!                ***  ROUTINE dom_vvl_ssh  ***
223      !!                   
224      !! ** Purpose :  compute the ssha field used in tra_zdf, dynnxt, tranxt
225      !!               and dynspg routines
226      !!----------------------------------------------------------------------
227      !! * Arguments
228      INTEGER, INTENT( in )  ::    kt                         ! time step
229
230      !! * Local declarations
231      INTEGER  ::   ji, jj, jk                                ! dummy loop indices
232      REAL(wp) ::   z2dt, zraur                               ! temporary scalars
233      REAL(wp), DIMENSION(jpi,jpj) ::   zun, zvn, zhdiv       ! 2D workspace
234      !!----------------------------------------------------------------------
235
236      !! Sea surface height after (ssha) in variable volume case
237      !! --------------------------====-----===============-----
238      !! ssha is needed for the tracer time-stepping (trazdf_imp or
239      !! tra_nxt), for the ssh time-stepping (dynspg_*) and for the
240      !! dynamics time-stepping (dynspg_flt or dyn_nxt, and wzv).
241      !! un, vn (or un_b and vn_b) and emp are needed, so it must be
242      !! done after the call of oce_sbc.
243
244      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
245         r2dt(:) =  rdttra(:)                       ! = rdtra (restarting with Euler time stepping)
246      ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1
247         r2dt(:) = 2. * rdttra(:)                   ! = 2 rdttra (leapfrog)
248      ENDIF
249
250      z2dt  = r2dt(1)
251      zraur = 1. / rauw
252
253      ! Vertically integrated quantities
254      ! --------------------------------
255      IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN
256         zun(:,:) = 0.e0    ;    zvn(:,:) = 0.e0
257         !
258         DO jk = 1, jpkm1
259            !                                               ! Vertically integrated transports (now)
260            zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk)
261            zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
262         END DO
263      ELSE ! lk_dynspg_ts is true
264         zun (:,:) = un_b(:,:)    ;    zvn (:,:) = vn_b(:,:)
265      ENDIF
266
267      ! Horizontal divergence of barotropic transports
268      !--------------------------------------------------
269      DO jj = 2, jpjm1
270         DO ji = 2, jpim1   ! vector opt.
271            zhdiv(ji,jj) = (  e2u(ji  ,jj  ) * zun(ji  ,jj  )    &
272               &            - e2u(ji-1,jj  ) * zun(ji-1,jj  )    &
273               &            + e1v(ji  ,jj  ) * zvn(ji  ,jj  )    &
274               &            - e1v(ji  ,jj-1) * zvn(ji  ,jj-1) )  &
275               &           / ( e1t(ji,jj) * e2t(ji,jj) )
276         END DO
277      END DO
278
279#if defined key_obc && ( key_dynspg_exp || key_dynspg_ts )
280      ! open boundaries (div must be zero behind the open boundary)
281      !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
282      IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east
283      IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west
284      IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
285      IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
286#endif
287
288      CALL lbc_lnk( zhdiv, 'T', 1. )
289
290      ! Sea surface elevation time stepping
291      ! -----------------------------------
292      IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN
293         ssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1)
294      ELSE ! lk_dynspg_ts is true
295         ssha(:,:) = (  sshb_b(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
296      ENDIF
297
298
299   END SUBROUTINE dom_vvl_ssh
300
301#else
302   !!----------------------------------------------------------------------
303   !!   Default option :                                      Empty routine
304   !!----------------------------------------------------------------------
305   SUBROUTINE dom_vvl_ini
306   END SUBROUTINE dom_vvl_ini
307   SUBROUTINE dom_vvl
308   END SUBROUTINE dom_vvl
309   SUBROUTINE dom_vvl_ssh( kt ) 
310      !! * Arguments
311      INTEGER, INTENT( in )  ::    kt        ! time step
312      WRITE(*,*) 'dom_vvl_ssh: You should not have seen this print! error?', kt
313   END SUBROUTINE dom_vvl_ssh
314#endif
315
316   !!======================================================================
317END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.