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 branches/2011/dev_r2802_NOCS_vvlfix/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2011/dev_r2802_NOCS_vvlfix/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 2827

Last change on this file since 2827 was 2827, checked in by acc, 13 years ago

Branch: dev_r2802_NOCS_vvlfix, minor correction to allow domvvl.F90 to be compiled without key_vvl

  • Property svn:keywords set to Id
File size: 11.8 KB
Line 
1MODULE domvvl
2   !!======================================================================
3   !!                       ***  MODULE domvvl   ***
4   !! Ocean :
5   !!======================================================================
6   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
7   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
8   !!----------------------------------------------------------------------
9#if defined key_vvl
10   !!----------------------------------------------------------------------
11   !!   'key_vvl'                              variable volume
12   !!----------------------------------------------------------------------
13   !!   dom_vvl     : defined coefficients to distribute ssh on each layers
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE sbc_oce         ! surface boundary condition: ocean
18   USE phycst          ! physical constants
19   USE in_out_manager  ! I/O manager
20   USE lib_mpp         ! distributed memory computing library
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   dom_vvl         ! called by domain.F90
27   PUBLIC   dom_vvl_2       ! called by domain.F90
28   PUBLIC   dom_vvl_alloc   ! called by nemogcm.F90
29
30   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mut , muu , muv , muf    !: 1/H_0 at t-,u-,v-,f-points
31
32   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:)     ::   r2dt   ! vertical profile time-step, = 2 rdttra
33      !                                                              ! except at nit000 (=rdttra) if neuler=0
34
35   !! * Substitutions
36#  include "domzgr_substitute.h90"
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS       
44
45   INTEGER FUNCTION dom_vvl_alloc()
46      !!----------------------------------------------------------------------
47      !!                ***  ROUTINE dom_vvl_alloc  ***
48      !!----------------------------------------------------------------------
49      !
50      ALLOCATE( mut (jpi,jpj,jpk) , muu (jpi,jpj,jpk) , muv (jpi,jpj,jpk) , muf (jpi,jpj,jpk) ,     &
51         &      r2dt        (jpk)                                                             , STAT=dom_vvl_alloc )
52         !
53      IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
54      IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
55      !
56   END FUNCTION dom_vvl_alloc
57
58
59   SUBROUTINE dom_vvl
60      !!----------------------------------------------------------------------
61      !!                ***  ROUTINE dom_vvl  ***
62      !!                   
63      !! ** Purpose :   compute mu coefficients at t-, u-, v- and f-points to
64      !!              spread ssh over the whole water column (scale factors)
65      !!                set the before and now ssh at u- and v-points
66      !!              (also f-point in now case)
67      !!----------------------------------------------------------------------
68      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
69      USE wrk_nemo, ONLY:   zee_t => wrk_2d_1, zee_u => wrk_2d_2, zee_v => wrk_2d_3, zee_f => wrk_2d_4   ! 2D workspace
70      !
71      INTEGER  ::   ji, jj, jk   ! dummy loop indices
72      REAL(wp) ::   zcoefu, zcoefv , zcoeff                ! local scalars
73      REAL(wp) ::   zvt   , zvt_ip1, zvt_jp1, zvt_ip1jp1   !   -      -
74      !!----------------------------------------------------------------------
75
76      IF( wrk_in_use(2, 1,2,3,4) ) THEN
77         CALL ctl_stop('dom_vvl: requested workspace arrays unavailable')   ;   RETURN
78      ENDIF
79
80      IF(lwp) THEN
81         WRITE(numout,*)
82         WRITE(numout,*) 'dom_vvl : Variable volume initialization'
83         WRITE(numout,*) '~~~~~~~~  compute coef. used to spread ssh over each layers'
84      ENDIF
85     
86      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl : unable to allocate arrays' )
87
88      fsdept(:,:,:) = gdept (:,:,:)
89      fsdepw(:,:,:) = gdepw (:,:,:)
90      fsde3w(:,:,:) = gdep3w(:,:,:)
91      fse3t (:,:,:) = e3t   (:,:,:)
92      fse3u (:,:,:) = e3u   (:,:,:)
93      fse3v (:,:,:) = e3v   (:,:,:)
94      fse3f (:,:,:) = e3f   (:,:,:)
95      fse3w (:,:,:) = e3w   (:,:,:)
96      fse3uw(:,:,:) = e3uw  (:,:,:)
97      fse3vw(:,:,:) = e3vw  (:,:,:)
98
99      !                                 !==  mu computation  ==!
100      zee_t(:,:) = fse3t_0(:,:,1)                ! Lower bound : thickness of the first model level
101      zee_u(:,:) = fse3u_0(:,:,1)
102      zee_v(:,:) = fse3v_0(:,:,1)
103      zee_f(:,:) = fse3f_0(:,:,1)
104      DO jk = 2, jpkm1                          ! Sum of the masked vertical scale factors
105         zee_t(:,:) = zee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk)
106         zee_u(:,:) = zee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
107         zee_v(:,:) = zee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
108         DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask
109            zee_f(:,jj) = zee_f(:,jj) + fse3f_0(:,jj,jk) *  umask(:,jj,jk) * umask(:,jj+1,jk)
110         END DO
111      END DO 
112      !                                         ! Compute and mask the inverse of the local depth at T, U, V and F points
113      zee_t(:,:) = 1._wp / zee_t(:,:) * tmask(:,:,1)
114      zee_u(:,:) = 1._wp / zee_u(:,:) * umask(:,:,1)
115      zee_v(:,:) = 1._wp / zee_v(:,:) * vmask(:,:,1)
116      DO jj = 1, jpjm1                               ! f-point case fmask cannot be used
117         zee_f(:,jj) = 1._wp / zee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1)
118      END DO
119      CALL lbc_lnk( zee_f, 'F', 1. )                 ! lateral boundary condition on ee_f
120      !
121      DO jk = 1, jpk                            ! mu coefficients
122         mut(:,:,jk) = zee_t(:,:) * tmask(:,:,jk)     ! T-point at T levels
123         muu(:,:,jk) = zee_u(:,:) * umask(:,:,jk)     ! U-point at T levels
124         muv(:,:,jk) = zee_v(:,:) * vmask(:,:,jk)     ! V-point at T levels
125      END DO
126      DO jk = 1, jpk                                 ! F-point : fmask=shlat at coasts, use the product of umask
127         DO jj = 1, jpjm1
128               muf(:,jj,jk) = zee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk)   ! at T levels
129         END DO
130         muf(:,jpj,jk) = 0._wp
131      END DO
132      CALL lbc_lnk( muf, 'F', 1. )                   ! lateral boundary condition
133
134
135      hu_0(:,:) = 0.e0                          ! Reference ocean depth at U- and V-points
136      hv_0(:,:) = 0.e0
137      DO jk = 1, jpk
138         hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
139         hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
140      END DO
141     
142      DO jj = 1, jpjm1                          ! initialise before and now Sea Surface Height at u-, v-, f-points
143         DO ji = 1, jpim1   ! NO vector opt.
144            zcoefu = 0.50_wp / ( e1u(ji,jj) * e2u(ji,jj) ) * umask(ji,jj,1)
145            zcoefv = 0.50_wp / ( e1v(ji,jj) * e2v(ji,jj) ) * vmask(ji,jj,1)
146            zcoeff = 0.25_wp / ( e1f(ji,jj) * e2f(ji,jj) ) * umask(ji,jj,1) * umask(ji,jj+1,1)
147            !
148            zvt           = e1e2t(ji  ,jj  ) * sshb(ji  ,jj  )    ! before fields
149            zvt_ip1       = e1e2t(ji+1,jj  ) * sshb(ji+1,jj  )
150            zvt_jp1       = e1e2t(ji  ,jj+1) * sshb(ji  ,jj+1)
151            sshu_b(ji,jj) = zcoefu * ( zvt + zvt_ip1 )
152            sshv_b(ji,jj) = zcoefv * ( zvt + zvt_jp1 )
153            !
154            zvt           = e1e2t(ji  ,jj  ) * sshn(ji  ,jj  )    ! now fields
155            zvt_ip1       = e1e2t(ji+1,jj  ) * sshn(ji+1,jj  )
156            zvt_jp1       = e1e2t(ji  ,jj+1) * sshn(ji  ,jj+1)
157            zvt_ip1jp1    = e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1)
158            sshu_n(ji,jj) = zcoefu * ( zvt + zvt_ip1 )
159            sshv_n(ji,jj) = zcoefv * ( zvt + zvt_jp1 )
160            sshf_n(ji,jj) = zcoeff * ( zvt + zvt_ip1 + zvt_jp1 + zvt_ip1jp1 )
161         END DO
162      END DO
163      CALL lbc_lnk( sshu_n, 'U', 1. )   ;   CALL lbc_lnk( sshu_b, 'U', 1. )      ! lateral boundary conditions
164      CALL lbc_lnk( sshv_n, 'V', 1. )   ;   CALL lbc_lnk( sshv_b, 'V', 1. )
165      CALL lbc_lnk( sshf_n, 'F', 1. )
166      !
167      IF( wrk_not_released(2, 1,2,3,4) )   CALL ctl_stop('dom_vvl: failed to release workspace arrays')
168      !
169   END SUBROUTINE dom_vvl
170
171
172   SUBROUTINE dom_vvl_2( kt, pe3u_b, pe3v_b )
173      !!----------------------------------------------------------------------
174      !!                ***  ROUTINE dom_vvl_2  ***
175      !!                   
176      !! ** Purpose :   compute the vertical scale factors at u- and v-points
177      !!              in variable volume case.
178      !!
179      !! ** Method  :   In variable volume case (non linear sea surface) the
180      !!              the vertical scale factor at velocity points is computed
181      !!              as the average of the cell surface weighted e3t.
182      !!                It uses the sea surface heigth so it have to be initialized
183      !!              after ssh is read/set
184      !!----------------------------------------------------------------------
185      INTEGER                   , INTENT(in   ) ::   kt               ! ocean time-step index
186      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe3u_b, pe3v_b   ! before vertical scale factor at u- & v-pts
187      !
188      INTEGER  ::   ji, jj, jk   ! dummy loop indices
189      INTEGER  ::   iku, ikv     ! local integers   
190      REAL(wp) ::   zvt          ! local scalars
191      !!----------------------------------------------------------------------
192
193      IF( lwp .AND. kt == nit000 ) THEN
194         WRITE(numout,*)
195         WRITE(numout,*) 'dom_vvl_2 : Variable volume, fse3t_b initialization'
196         WRITE(numout,*) '~~~~~~~~~ '
197         pe3u_b(:,:,jpk) = fse3u_0(:,:,jpk)
198         pe3v_b(:,:,jpk) = fse3u_0(:,:,jpk)
199      ENDIF
200     
201      DO jk = 1, jpkm1           ! set the before scale factors at u- & v-points
202         DO jj = 2, jpjm1
203            DO ji = fs_2, fs_jpim1
204               zvt = fse3t_b(ji,jj,jk) * e1e2t(ji,jj)
205               pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1e2t(ji+1,jj) ) / ( e1u(ji,jj) * e2u(ji,jj) )
206               pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e1e2t(ji,jj+1) ) / ( e1v(ji,jj) * e2v(ji,jj) )
207            END DO
208         END DO
209      END DO
210      IF( ln_zps ) THEN          ! minimum of the e3t at partial cell level
211         DO jj = 2, jpjm1
212            DO ji = fs_2, fs_jpim1
213               iku = mbku(ji,jj)
214               ikv = mbkv(ji,jj)
215               pe3u_b(ji,jj,iku) = MIN( fse3t_b(ji,jj,iku), fse3t_b(ji+1,jj  ,iku) ) 
216               pe3v_b(ji,jj,ikv) = MIN( fse3t_b(ji,jj,ikv), fse3t_b(ji  ,jj+1,ikv) ) 
217            END DO
218         END DO
219      ENDIF
220      pe3u_b(:,:,:) = pe3u_b(:,:,:) - fse3u_0(:,:,:)      ! anomaly to avoid zero along closed boundary/extra halos
221      pe3v_b(:,:,:) = pe3v_b(:,:,:) - fse3v_0(:,:,:)
222      CALL lbc_lnk( pe3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions
223      CALL lbc_lnk( pe3v_b(:,:,:), 'V', 1. )
224      pe3u_b(:,:,:) = pe3u_b(:,:,:) + fse3u_0(:,:,:)      ! recover the full scale factor
225      pe3v_b(:,:,:) = pe3v_b(:,:,:) + fse3v_0(:,:,:)
226      !
227   END SUBROUTINE dom_vvl_2
228   
229#else
230   !!----------------------------------------------------------------------
231   !!   Default option :                                      Empty routine
232   !!----------------------------------------------------------------------
233CONTAINS
234   SUBROUTINE dom_vvl
235   END SUBROUTINE dom_vvl
236   SUBROUTINE dom_vvl_2(kdum, pudum, pvdum )
237      USE par_kind
238      INTEGER                   , INTENT(in   ) ::   kdum       
239      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pudum, pvdum
240   END SUBROUTINE dom_vvl_2
241#endif
242
243   !!======================================================================
244END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.