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_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 3432

Last change on this file since 3432 was 3432, checked in by trackstand2, 12 years ago

Merge branch 'ksection_partition'
  • Property svn:keywords set to Id
File size: 10.6 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_alloc ! called by nemogcm.F90
28
29   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ee_t, ee_u, ee_v, ee_f   !: ???
30   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mut , muu , muv , muf    !: ???
31
32   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:)     ::   r2dt   ! vertical profile time-step, = 2 rdttra
33      !                                                              ! except at nit000 (=rdttra) if neuler=0
34
35   !! * Control permutation of array indices
36#  include "oce_ftrans.h90"
37#  include "dom_oce_ftrans.h90"
38#  include "sbc_oce_ftrans.h90"
39#  include "domvvl_ftrans.h90"
40
41   !! * Substitutions
42#  include "domzgr_substitute.h90"
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
46   !! $Id$
47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49CONTAINS       
50
51   INTEGER FUNCTION dom_vvl_alloc()
52      !!----------------------------------------------------------------------
53      !!                ***  ROUTINE dom_vvl_alloc  ***
54      !!----------------------------------------------------------------------
55      !
56      ALLOCATE( mut (jpi,jpj,jpk) , muu (jpi,jpj,jpk) , muv (jpi,jpj,jpk) , muf (jpi,jpj,jpk) ,     &
57         &      ee_t(jpi,jpj)     , ee_u(jpi,jpj)     , ee_v(jpi,jpj)     , ee_f(jpi,jpj)     ,     &
58         &      r2dt        (jpk)                                                             , STAT=dom_vvl_alloc )
59         !
60      IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
61      IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
62      !
63   END FUNCTION dom_vvl_alloc
64
65
66   SUBROUTINE dom_vvl
67      !!----------------------------------------------------------------------
68      !!                ***  ROUTINE dom_vvl  ***
69      !!                   
70      !! ** Purpose :  compute coefficients muX at T-U-V-F points to spread
71      !!               ssh over the whole water column (scale factors)
72      !!----------------------------------------------------------------------
73      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
74      USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3     ! 2D workspace
75      !
76      INTEGER  ::   ji, jj, jk   ! dummy loop indices
77      REAL(wp) ::   zcoefu , zcoefv   , zcoeff                   ! local scalars
78      REAL(wp) ::   zv_t_ij, zv_t_ip1j, zv_t_ijp1, zv_t_ip1jp1   !   -      -
79      !!----------------------------------------------------------------------
80
81      IF( wrk_in_use(2, 1,2,3) ) THEN
82         CALL ctl_stop('dom_vvl: requested workspace arrays unavailable')   ;   RETURN
83      ENDIF
84
85      IF(lwp) THEN
86         WRITE(numout,*)
87         WRITE(numout,*) 'dom_vvl : Variable volume initialization'
88         WRITE(numout,*) '~~~~~~~~  compute coef. used to spread ssh over each layers'
89      ENDIF
90     
91      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl : unable to allocate arrays' )
92
93      fsdept(:,:,:) = gdept (:,:,:)
94      fsdepw(:,:,:) = gdepw (:,:,:)
95      fsde3w(:,:,:) = gdep3w(:,:,:)
96      fse3t (:,:,:) = e3t   (:,:,:)
97      fse3u (:,:,:) = e3u   (:,:,:)
98      fse3v (:,:,:) = e3v   (:,:,:)
99      fse3f (:,:,:) = e3f   (:,:,:)
100      fse3w (:,:,:) = e3w   (:,:,:)
101      fse3uw(:,:,:) = e3uw  (:,:,:)
102      fse3vw(:,:,:) = e3vw  (:,:,:)
103
104      !                                 !==  mu computation  ==!
105      ee_t(:,:) = fse3t_0(:,:,1)                ! Lower bound : thickness of the first model level
106      ee_u(:,:) = fse3u_0(:,:,1)
107      ee_v(:,:) = fse3v_0(:,:,1)
108      ee_f(:,:) = fse3f_0(:,:,1)
109      DO jk = 2, jpkm1                          ! Sum of the masked vertical scale factors
110         ee_t(:,:) = ee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk)
111         ee_u(:,:) = ee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
112         ee_v(:,:) = ee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
113         DO jj = 1, jpjm1                      ! f-point : fmask=shlat at coasts, use the product of umask
114            ee_f(:,jj) = ee_f(:,jj) + fse3f_0(:,jj,jk) *  umask(:,jj,jk) * umask(:,jj+1,jk)
115         END DO
116      END DO 
117      !                                         ! Compute and mask the inverse of the local depth at T, U, V and F points
118#if defined key_z_first
119      ee_t(:,:) = 1. / ee_t(:,:) * tmask_1(:,:)
120      ee_u(:,:) = 1. / ee_u(:,:) * umask_1(:,:)
121      ee_v(:,:) = 1. / ee_v(:,:) * vmask_1(:,:)
122      DO jj = 1, jpjm1                               ! f-point case fmask cannot be used
123         ee_f(:,jj) = 1. / ee_f(:,jj) * umask_1(:,jj) * umask_1(:,jj+1)
124      END DO
125#else
126      ee_t(:,:) = 1. / ee_t(:,:) * tmask(:,:,1)
127      ee_u(:,:) = 1. / ee_u(:,:) * umask(:,:,1)
128      ee_v(:,:) = 1. / ee_v(:,:) * vmask(:,:,1)
129      DO jj = 1, jpjm1                               ! f-point case fmask cannot be used
130         ee_f(:,jj) = 1. / ee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1)
131      END DO
132#endif
133      CALL lbc_lnk( ee_f, 'F', 1. )                  ! lateral boundary condition on ee_f
134      !
135      DO jk = 1, jpk                            ! mu coefficients
136         mut(:,:,jk) = ee_t(:,:) * tmask(:,:,jk)     ! T-point at T levels
137         muu(:,:,jk) = ee_u(:,:) * umask(:,:,jk)     ! U-point at T levels
138         muv(:,:,jk) = ee_v(:,:) * vmask(:,:,jk)     ! V-point at T levels
139      END DO
140      DO jk = 1, jpk                                 ! F-point : fmask=shlat at coasts, use the product of umask
141         DO jj = 1, jpjm1
142               muf(:,jj,jk) = ee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk)   ! at T levels
143         END DO
144         muf(:,jpj,jk) = 0.e0
145      END DO
146      CALL lbc_lnk( muf, 'F', 1. )                   ! lateral boundary condition
147
148
149      hu_0(:,:) = 0.e0                          ! Reference ocean depth at U- and V-points
150      hv_0(:,:) = 0.e0
151      DO jk = 1, jpk
152         hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
153         hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
154      END DO
155     
156      ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
157      ! for ssh and scale factors
158      zs_t  (:,:) =         e1t(:,:) * e2t(:,:)
159      zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )
160      zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )
161
162      DO jj = 1, jpjm1                          ! initialise before and now Sea Surface Height at u-, v-, f-points
163         DO ji = 1, jpim1   ! NO vector opt.
164            zcoefu = umask(ji,jj,1) * zs_u_1(ji,jj)
165            zcoefv = vmask(ji,jj,1) * zs_v_1(ji,jj)
166            zcoeff = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) / ( e1f(ji,jj) * e2f(ji,jj) )
167            ! before fields
168            zv_t_ij       = zs_t(ji  ,jj  ) * sshb(ji  ,jj  )
169            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshb(ji+1,jj  )
170            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshb(ji  ,jj+1)
171            sshu_b(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j )
172            sshv_b(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 )
173            ! now fields
174            zv_t_ij       = zs_t(ji  ,jj  ) * sshn(ji  ,jj  )
175            zv_t_ip1j     = zs_t(ji+1,jj  ) * sshn(ji+1,jj  )
176            zv_t_ijp1     = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1)
177            zv_t_ip1jp1   = zs_t(ji  ,jj+1) * sshn(ji  ,jj+1)
178            sshu_n(ji,jj) = zcoefu * ( zv_t_ij + zv_t_ip1j )
179            sshv_n(ji,jj) = zcoefv * ( zv_t_ij + zv_t_ijp1 )
180            sshf_n(ji,jj) = zcoeff * ( zv_t_ij + zv_t_ip1j + zv_t_ijp1 + zv_t_ip1jp1 )
181         END DO
182      END DO
183      CALL lbc_lnk( sshu_n, 'U', 1. )   ;   CALL lbc_lnk( sshu_b, 'U', 1. )      ! lateral boundary conditions
184      CALL lbc_lnk( sshv_n, 'V', 1. )   ;   CALL lbc_lnk( sshv_b, 'V', 1. )
185      CALL lbc_lnk( sshf_n, 'F', 1. )
186
187                                                ! initialise before scale factors at (u/v)-points
188      ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
189#if defined key_z_first
190      DO jj = 1, jpjm1
191         DO ji = 1, jpim1
192            DO jk = 1, jpkm1
193#else
194      DO jk = 1, jpkm1
195         DO jj = 1, jpjm1
196            DO ji = 1, jpim1
197#endif
198               zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
199               zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
200               zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
201               fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
202               fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
203            END DO
204         END DO
205      END DO
206      CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )               ! lateral boundary conditions
207      CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
208      ! Add initial scale factor to scale factor anomaly
209      fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
210      fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
211      !
212      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dom_vvl: failed to release workspace arrays')
213      !
214   END SUBROUTINE dom_vvl
215
216#else
217   !!----------------------------------------------------------------------
218   !!   Default option :                                      Empty routine
219   !!----------------------------------------------------------------------
220CONTAINS
221   SUBROUTINE dom_vvl
222   END SUBROUTINE dom_vvl
223#endif
224
225   !!======================================================================
226END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.