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

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • 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.