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

Last change on this file since 1058 was 1058, checked in by rblod, 16 years ago

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