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 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.1 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   PUBLIC dom_vvl_sf_ini !
35   PUBLIC dom_vvl_sf     !
36
37   !! * Module variables
38   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !:
39      mut, muu, muv, muf                            !:
40
41   REAL(wp), DIMENSION(jpk) ::   r2dt               ! vertical profile time-step, = 2 rdttra
42      !                                             ! except at nit000 (=rdttra) if neuler=0
43
44   !! * Substitutions
45#  include "domzgr_substitute.h90"
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !!   OPA 9.0 , LOCEAN-IPSL (2005)
49   !! $Header$
50   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52
53CONTAINS       
54
55#if defined key_vvl
56
57   SUBROUTINE dom_vvl_ini
58      !!----------------------------------------------------------------------
59      !!                ***  ROUTINE dom_vvl_ini  ***
60      !!                   
61      !! ** Purpose :  compute coefficients muX at T-U-V-F points to spread
62      !!               ssh over the whole water column (scale factors)
63      !!
64      !!----------------------------------------------------------------------
65      INTEGER  ::   ji, jj, jk, zpk
66      REAL(wp), DIMENSION(jpi,jpj) ::   zmut, zmuu, zmuv, zmuf   ! 2D workspace
67      !!----------------------------------------------------------------------
68
69      IF(lwp)   THEN
70         WRITE(numout,*)
71         WRITE(numout,*) 'dom_vvl_ini : Variable volume activated'
72         WRITE(numout,*) '~~~~~~~~~~~   compute coef. used to spread ssh over each layers'
73      ENDIF
74
75      IF( ln_zps )  CALL ctl_stop( 'dom_vvl_ini : option ln_zps is incompatible with variable volume option key_vvl')
76
77#if defined key_zco  ||  defined key_dynspg_rl
78      CALL ctl_stop( 'dom_vvl_ini : options key_zco/key_dynspg_rl are incompatible with variable volume option key_vvl')
79#endif
80
81      fsvdept (:,:,:) = gdept (:,:,:)
82      fsvdepw (:,:,:) = gdepw (:,:,:)
83      fsvde3w (:,:,:) = gdep3w(:,:,:)
84      fsve3t  (:,:,:) = e3t   (:,:,:)
85      fsve3u  (:,:,:) = e3u   (:,:,:)
86      fsve3v  (:,:,:) = e3v   (:,:,:)
87      fsve3f  (:,:,:) = e3f   (:,:,:)
88      fsve3w  (:,:,:) = e3w   (:,:,:)
89      fsve3uw (:,:,:) = e3uw  (:,:,:)
90      fsve3vw (:,:,:) = e3vw  (:,:,:)
91
92      ! mu computation
93      zmut(:,:)   = 0.e0
94      zmuu(:,:)   = 0.e0
95      zmuv(:,:)   = 0.e0
96      zmuf(:,:)   = 0.e0
97      mut (:,:,:) = 0.e0
98      muu (:,:,:) = 0.e0
99      muv (:,:,:) = 0.e0
100      muf (:,:,:) = 0.e0
101
102      DO jj = 1, jpj
103         DO ji = 1, jpi
104            zpk = mbathy(ji,jj) - 1
105            DO jk = 1, zpk
106               zmut(ji,jj) = zmut(ji,jj) + fsve3t(ji,jj,jk) * SUM( fsve3t(ji,jj,jk:zpk) )
107               zmuu(ji,jj) = zmuu(ji,jj) + fsve3u(ji,jj,jk) * SUM( fsve3u(ji,jj,jk:zpk) )
108               zmuv(ji,jj) = zmuv(ji,jj) + fsve3v(ji,jj,jk) * SUM( fsve3v(ji,jj,jk:zpk) )
109               zmuf(ji,jj) = zmuf(ji,jj) + fsve3f(ji,jj,jk) * SUM( fsve3f(ji,jj,jk:zpk) )
110            END DO
111         END DO
112      END DO
113
114      DO jj = 1, jpj
115         DO ji = 1, jpi
116            zpk = mbathy(ji,jj) - 1
117            DO jk = 1, zpk
118               mut(ji,jj,jk) = SUM( fsve3t(ji,jj,jk:zpk) ) / zmut(ji,jj)
119               muu(ji,jj,jk) = SUM( fsve3u(ji,jj,jk:zpk) ) / zmuu(ji,jj)
120               muv(ji,jj,jk) = SUM( fsve3v(ji,jj,jk:zpk) ) / zmuv(ji,jj)
121               muf(ji,jj,jk) = SUM( fsve3f(ji,jj,jk:zpk) ) / zmuf(ji,jj)
122            END DO
123         END DO
124      END DO
125
126
127   END SUBROUTINE dom_vvl_ini
128
129
130
131   SUBROUTINE dom_vvl
132      !!----------------------------------------------------------------------
133      !!                ***  ROUTINE dom_vvl  ***
134      !!                   
135      !! ** Purpose :  compute ssh at U-V-F points, T-W scale factors and local
136      !!               depths at each time step.
137      !!----------------------------------------------------------------------
138      !! * Local declarations
139      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
140      REAL(wp), DIMENSION(jpi,jpj) :: zsshf    ! 2D workspace
141      !!----------------------------------------------------------------------
142
143      ! Sea Surface Height at u-v-fpoints
144      DO jj = 1, jpjm1
145         DO ji = 1,jpim1
146            sshu(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )  &
147               &          * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
148               &            + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
149            !
150            sshv(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )  &
151               &          * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
152               &            + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
153            !
154            zsshf(ji,jj) = 0.25 * fmask(ji,jj,1)                 &
155               &           * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   &
156               &             + sshn(ji+1,jj) + sshn(ji+1,jj+1) )
157         END DO
158      END DO
159
160      ! Boundaries conditions
161      CALL lbc_lnk( sshu,  'U', 1. )
162      CALL lbc_lnk( sshv,  'V', 1. )
163      CALL lbc_lnk( zsshf, 'F', 1. )
164
165      ! Scale factors at T levels
166      DO jk = 1, jpkm1
167         fse3t(:,:,jk) = fsve3t(:,:,jk) * ( 1 + sshn(:,:)  * mut(:,:,jk) )
168         fse3u(:,:,jk) = fsve3u(:,:,jk) * ( 1 + sshu(:,:)  * muu(:,:,jk) )
169         fse3v(:,:,jk) = fsve3v(:,:,jk) * ( 1 + sshv(:,:)  * muv(:,:,jk) )
170         fse3f(:,:,jk) = fsve3f(:,:,jk) * ( 1 + zsshf(:,:) * muf(:,:,jk) )
171      END DO
172
173      ! Scale factors at W levels
174      fse3w (:,:,1) = fse3t(:,:,1)
175      fse3uw(:,:,1) = fse3u(:,:,1)
176      fse3vw(:,:,1) = fse3v(:,:,1)
177      DO jk = 2, jpk
178         fse3w (:,:,jk) = 0.5 * ( fse3t(:,:,jk-1) + fse3t(:,:,jk) )
179         fse3uw(:,:,jk) = 0.5 * ( fse3u(:,:,jk-1) + fse3u(:,:,jk) )
180         fse3vw(:,:,jk) = 0.5 * ( fse3v(:,:,jk-1) + fse3v(:,:,jk) )
181      END DO
182
183      ! T and W points depth
184      fsdept(:,:,1) = 0.5 * fse3w(:,:,1)
185      fsdepw(:,:,1) = 0.e0
186      fsde3w(:,:,1) = fsdept(:,:,1) - sshn(:,:)
187      DO jk = 2, jpk
188         fsdept(:,:,jk) = fsdept(:,:,jk-1) + fse3w(:,:,jk)
189         fsdepw(:,:,jk) = fsdepw(:,:,jk-1) + fse3t(:,:,jk-1)
190         fsde3w(:,:,jk) = fsdept(:,:,jk  ) - sshn(:,:)
191      END DO
192
193      ! Local depth or Inverse of the local depth of the water column at u- and v-points
194      ! ------------------------------
195
196      ! Ocean depth at U- and V-points
197      hu(:,:) = 0.e0    ;    hv(:,:) = 0.e0
198
199      DO jk = 1, jpk
200         hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk)
201         hv(:,:) = hv(:,:) + fse3v(:,:,jk) * vmask(:,:,jk)
202      END DO
203
204
205      ! Inverse of the local depth
206      hur(:,:) = fse3u(:,:,1)             ! Lower bound : thickness of the first model level
207      hvr(:,:) = fse3v(:,:,1)
208     
209      DO jk = 2, jpk                      ! Sum of the vertical scale factors
210         hur(:,:) = hur(:,:) + fse3u(:,:,jk) * umask(:,:,jk)
211         hvr(:,:) = hvr(:,:) + fse3v(:,:,jk) * vmask(:,:,jk)
212      END DO
213
214      ! Compute and mask the inverse of the local depth
215      hur(:,:) = 1. / hur(:,:) * umask(:,:,1)
216      hvr(:,:) = 1. / hvr(:,:) * vmask(:,:,1)
217
218   END SUBROUTINE dom_vvl
219
220
221
222   SUBROUTINE dom_vvl_ssh( kt ) 
223      !!----------------------------------------------------------------------
224      !!                ***  ROUTINE dom_vvl_ssh  ***
225      !!                   
226      !! ** Purpose :  compute the ssha field used in tra_zdf, dynnxt, tranxt
227      !!               and dynspg routines
228      !!----------------------------------------------------------------------
229      !! * Arguments
230      INTEGER, INTENT( in )  ::    kt                         ! time step
231
232      !! * Local declarations
233      INTEGER  ::   ji, jj, jk                                ! dummy loop indices
234      REAL(wp) ::   z2dt, zraur                               ! temporary scalars
235      REAL(wp), DIMENSION(jpi,jpj) ::   zun, zvn, zhdiv       ! 2D workspace
236      !!----------------------------------------------------------------------
237
238      !! Sea surface height after (ssha) in variable volume case
239      !! --------------------------====-----===============-----
240      !! ssha is needed for the tracer time-stepping (trazdf_imp or
241      !! tra_nxt), for the ssh time-stepping (dynspg_*) and for the
242      !! dynamics time-stepping (dynspg_flt or dyn_nxt, and wzv).
243      !! un, vn (or un_b and vn_b) and emp are needed, so it must be
244      !! done after the call of oce_sbc.
245
246      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
247         r2dt(:) =  rdttra(:)                       ! = rdtra (restarting with Euler time stepping)
248      ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1
249         r2dt(:) = 2. * rdttra(:)                   ! = 2 rdttra (leapfrog)
250      ENDIF
251
252      z2dt  = r2dt(1)
253      zraur = 1. / rauw
254
255      ! Vertically integrated quantities
256      ! --------------------------------
257      IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN
258         zun(:,:) = 0.e0    ;    zvn(:,:) = 0.e0
259         !
260         DO jk = 1, jpkm1
261            !                                               ! Vertically integrated transports (now)
262            zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk)
263            zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
264         END DO
265      ELSE ! lk_dynspg_ts is true
266         zun (:,:) = un_b(:,:)    ;    zvn (:,:) = vn_b(:,:)
267      ENDIF
268
269      ! Horizontal divergence of barotropic transports
270      !--------------------------------------------------
271      DO jj = 2, jpjm1
272         DO ji = 2, jpim1   ! vector opt.
273            zhdiv(ji,jj) = (  e2u(ji  ,jj  ) * zun(ji  ,jj  )    &
274               &            - e2u(ji-1,jj  ) * zun(ji-1,jj  )    &
275               &            + e1v(ji  ,jj  ) * zvn(ji  ,jj  )    &
276               &            - e1v(ji  ,jj-1) * zvn(ji  ,jj-1) )  &
277               &           / ( e1t(ji,jj) * e2t(ji,jj) )
278         END DO
279      END DO
280
281#if defined key_obc && ( key_dynspg_exp || key_dynspg_ts )
282      ! open boundaries (div must be zero behind the open boundary)
283      !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
284      IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east
285      IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west
286      IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
287      IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
288#endif
289
290      CALL lbc_lnk( zhdiv, 'T', 1. )
291
292      ! Sea surface elevation time stepping
293      ! -----------------------------------
294      IF( .NOT. lk_dynspg_ts .OR. (neuler == 0 .AND. kt == nit000) ) THEN
295         ssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1)
296      ELSE ! lk_dynspg_ts is true
297         ssha(:,:) = (  sshb_b(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
298      ENDIF
299
300
301   END SUBROUTINE dom_vvl_ssh
302
303   SUBROUTINE  dom_vvl_sf( zssh, gridp, sfe3 )
304      !!----------------------------------------------------------------------
305      !!                ***  ROUTINE dom_vvl_sf  ***
306      !!                   
307      !! ** Purpose :  compute vertical scale factor at each time step
308      !!----------------------------------------------------------------------
309      !! * Arguments
310      CHARACTER(LEN=1)                , INTENT( in ) :: gridp   ! grid point type
311      REAL(wp), DIMENSION(jpi,jpj)    , INTENT( in ) :: zssh    ! 2D workspace
312      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: sfe3    ! 3D workspace
313
314      !! * Local declarations
315      INTEGER  ::   jk                                      ! dummy loop indices
316      !!----------------------------------------------------------------------
317
318      SELECT CASE (gridp)
319
320      CASE ('U')
321         !
322         DO jk = 1, jpk
323            sfe3(:,:,jk)  = fsve3u(:,:,jk) * ( 1 + zssh(:,:) * muu(:,:,jk) )
324         ENDDO
325
326      CASE ('V')
327         !
328         DO jk = 1, jpk
329            sfe3(:,:,jk)  = fsve3v(:,:,jk) * ( 1 + zssh(:,:) * muv(:,:,jk) )
330         ENDDO
331
332      CASE ('T')
333         !
334         DO jk = 1, jpk
335            sfe3(:,:,jk)  = fsve3t(:,:,jk) * ( 1 + zssh(:,:) * mut(:,:,jk) )
336         ENDDO
337
338      END SELECT
339
340   END SUBROUTINE dom_vvl_sf
341
342   SUBROUTINE dom_vvl_sf_ini( gridp, sfe3ini )
343      !!----------------------------------------------------------------------
344      !!                ***  ROUTINE dom_vvl_sf_ini  ***
345      !!                   
346      !! ** Purpose :  affect the appropriate vertical scale factor. It is done
347      !!               to avoid compiling error when using key_zco cpp key
348      !!----------------------------------------------------------------------
349      !! * Arguments
350      CHARACTER(LEN=1)                , INTENT( in ) :: gridp      ! grid point type
351      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3ini    ! 3D workspace
352      !!----------------------------------------------------------------------
353
354      SELECT CASE (gridp)
355
356      CASE ('U')
357         !
358         sfe3ini(:,:,:)  = fse3u(:,:,:)
359
360      CASE ('V')
361         !
362         sfe3ini(:,:,:)  = fse3v(:,:,:)
363
364      CASE ('T')
365         !
366         sfe3ini(:,:,:)  = fse3t(:,:,:)
367
368      END SELECT
369
370   END SUBROUTINE dom_vvl_sf_ini
371#else
372   !!----------------------------------------------------------------------
373   !!   Default option :                                      Empty routine
374   !!----------------------------------------------------------------------
375   SUBROUTINE dom_vvl_ini
376   END SUBROUTINE dom_vvl_ini
377   SUBROUTINE dom_vvl
378   END SUBROUTINE dom_vvl
379   SUBROUTINE dom_vvl_ssh( kt ) 
380      !! * Arguments
381      INTEGER, INTENT( in )  ::    kt        ! time step
382      WRITE(*,*) 'dom_vvl_ssh: You should not have seen this print! error?', kt
383   END SUBROUTINE dom_vvl_ssh
384   SUBROUTINE dom_vvl_sf( zssh, gridp, sfe3 ) 
385      !! * Arguments
386      CHARACTER(LEN=1)                , INTENT( in ) :: gridp   ! grid point type
387      REAL(wp), DIMENSION(jpi,jpj)    , INTENT( in ) :: zssh    ! 2D workspace
388      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3    ! 3D workspace
389      sfe3(:,:,:) = 0.e0
390      WRITE(*,*) 'sfe3: You should not have seen this print! error?', gridp
391      WRITE(*,*) 'sfe3: You should not have seen this print! error?', zssh(1,1)
392   END SUBROUTINE dom_vvl_sf
393   SUBROUTINE dom_vvl_sf_ini( gridp, sfe3ini ) 
394      !! * Arguments
395      CHARACTER(LEN=1)                , INTENT( in ) :: gridp    ! grid point type
396      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT (out) :: sfe3ini  ! 3D workspace
397      sfe3ini(:,:,:) = 0.e0
398      WRITE(*,*) 'sfe3ini: You should not have seen this print! error?', gridp
399   END SUBROUTINE dom_vvl_sf_ini
400#endif
401
402   !!======================================================================
403END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.