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.
domqco.F90 in NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/tests/VORTEX/MY_SRC – NEMO

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/tests/VORTEX/MY_SRC/domqco.F90 @ 13679

Last change on this file since 13679 was 13679, checked in by jchanut, 3 years ago

#2345. Add volume initialization with qco in VORTEX test case. At this stage, test case compiles and and runs but parent grid surface kinematic condition appears to be wrong under overlap.
File size: 25.5 KB
Line 
1MODULE domqco
2   !!======================================================================
3   !!                       ***  MODULE domqco   ***
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   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates
9   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
10   !!            4.1  !  2019-08  (A. Coward, D. Storkey) add time level indices for prognostic variables
11   !!            4.x  !  2020-02  (S. Techene, G. Madec) quasi-eulerian coordinate (z* or s*) from domvvl
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   dom_qco_init   : define initial vertical scale factors, depths and column thickness
16   !!   dom_qco_zgr    : Set ssh/h_0 ratio at t
17   !!   dom_qco_r3c    : Compute ssh/h_0 ratio at t-, u-, v-, and optionally f-points
18   !!       qco_rst_read : read/write restart file
19   !!       qco_ctl    : Check the vvl options
20   !!----------------------------------------------------------------------
21   USE oce            ! ocean dynamics and tracers
22   USE phycst         ! physical constant
23   USE dom_oce        ! ocean space and time domain
24   USE dynadv  , ONLY : ln_dynadv_vec
25   USE isf_oce        ! iceshelf cavities
26   USE sbc_oce        ! ocean surface boundary condition
27   USE wet_dry        ! wetting and drying
28   USE usrdef_istate  ! user defined initial state (wad only)
29   USE restart        ! ocean restart
30   !
31   USE in_out_manager ! I/O manager
32   USE iom            ! I/O manager library
33   USE lib_mpp        ! distributed memory computing library
34   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
35   USE timing         ! Timing
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC  dom_qco_init       ! called by domain.F90
41   PUBLIC  dom_qco_zgr        ! called by isfcpl.F90
42   PUBLIC  dom_qco_r3c        ! called by steplf.F90
43
44   !                                                      !!* Namelist nam_vvl
45   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate
46   LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate
47   LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate
48   LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
49   LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
50   LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
51   !                                                       ! conservation: not used yet
52   REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient
53   REAL(wp)         :: rn_rst_e3t                          ! ztilde to zstar restoration timescale [days]
54   REAL(wp)         :: rn_lf_cutoff                        ! cutoff frequency for low-pass filter  [days]
55   REAL(wp)         :: rn_zdef_max                         ! maximum fractional e3t deformation
56   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints
57
58   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport
59
60   !! * Substitutions
61#  include "do_loop_substitute.h90"
62   !!----------------------------------------------------------------------
63   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
64   !! $Id: domvvl.F90 12377 2020-02-12 14:39:06Z acc $
65   !! Software governed by the CeCILL license (see ./LICENSE)
66   !!----------------------------------------------------------------------
67CONTAINS
68
69   SUBROUTINE dom_qco_init( Kbb, Kmm, Kaa )
70      !!----------------------------------------------------------------------
71      !!                ***  ROUTINE dom_qco_init  ***
72      !!
73      !! ** Purpose :  Initialization of all ssh. to h._0 ratio
74      !!
75      !! ** Method  :  - use restart file and/or initialize
76      !!               - compute ssh. to h._0 ratio
77      !!
78      !! ** Action  : - r3(t/u/v)_b
79      !!              - r3(t/u/v/f)_n
80      !!
81      !!----------------------------------------------------------------------
82      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! time level indices
83      !!----------------------------------------------------------------------
84      !
85      IF(lwp) WRITE(numout,*)
86      IF(lwp) WRITE(numout,*) 'dom_qco_init : Variable volume activated'
87      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
88      !
89      CALL qco_ctl                            ! choose vertical coordinate (z_star, z_tilde or layer)
90      !              ! CAUTION COM A METTRE !!!
91!!st      CALL qco_rst_read2( nit000, Kbb, Kmm )  ! Read or initialize ssh_(Kbb/Kmm) and r3
92!!st CAUTION if read2 removed change restart.F90 !
93      !
94      CALL qco_rst_read( nit000, Kbb, Kmm )   ! Read or initialize ssh_(Kbb/Kmm)
95      !
96      CALL dom_qco_zgr( Kbb, Kmm, Kaa )       ! interpolation scale factor, depth and water column
97      !
98      IF(lwxios) THEN   ! define variables in restart file when writing with XIOS
99         CALL iom_set_rstw_var_active('sshb')
100         CALL iom_set_rstw_var_active('sshn')
101      ENDIF
102      !
103   END SUBROUTINE dom_qco_init
104
105
106   SUBROUTINE dom_qco_zgr( Kbb, Kmm, Kaa )
107      !!----------------------------------------------------------------------
108      !!                ***  ROUTINE dom_qco_init  ***
109      !!
110      !! ** Purpose :  Initialization of all ssh./h._0 ratio
111      !!
112      !! ** Method  :  - call domqco using Kbb and Kmm
113      !!
114      !! ** Action  : - r3(t/u/v)_b
115      !!              - r3(t/u/v/f)_n
116      !!----------------------------------------------------------------------
117      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! time level indices
118      !!----------------------------------------------------------------------
119      !
120      !                    !== Set of all other vertical scale factors  ==!  (now and before)
121      !                                ! Horizontal interpolation of e3t
122      CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb) )
123      CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) )
124      !
125#if defined key_agrif
126      ! We need to define r3[tuv](Kaa) for AGRIF initialisation (should not be a problem for the restartability...)
127      r3t(:,:,Kaa) = r3t(:,:,Kmm)
128      r3u(:,:,Kaa) = r3u(:,:,Kmm)
129      r3v(:,:,Kaa) = r3v(:,:,Kmm)
130#endif
131      !
132   END SUBROUTINE dom_qco_zgr
133
134
135   SUBROUTINE dom_qco_r3c( pssh, pr3t, pr3u, pr3v, pr3f )
136      !!---------------------------------------------------------------------
137      !!                   ***  ROUTINE r3c  ***
138      !!
139      !! ** Purpose :   compute the filtered ratio ssh/h_0 at t-,u-,v-,f-points
140      !!
141      !! ** Method  : - compute the ssh at u- and v-points (f-point optional)
142      !!                   Vector Form : surface weighted averaging
143      !!                   Flux   Form : simple           averaging
144      !!              - compute the ratio ssh/h_0 at t-,u-,v-pts, (f-pt optional)
145      !!----------------------------------------------------------------------
146      REAL(wp), DIMENSION(:,:)          , INTENT(in   )  ::   pssh               ! sea surface height   [m]
147      REAL(wp), DIMENSION(:,:)          , INTENT(  out)  ::   pr3t, pr3u, pr3v   ! ssh/h0 ratio at t-, u-, v-,points  [-]
148      REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT(  out)  ::   pr3f               ! ssh/h0 ratio at f-point   [-]
149      !
150      INTEGER ::   ji, jj   ! dummy loop indices
151      !!----------------------------------------------------------------------
152      !
153      !
154      pr3t(:,:) = pssh(:,:) * r1_ht_0(:,:)   !==  ratio at t-point  ==!
155      !
156      !
157      !                                      !==  ratio at u-,v-point  ==!
158      !
159!!st      IF( ln_dynadv_vec ) THEN                     !- Vector Form   (thickness weighted averaging)
160#if ! defined key_qcoTest_FluxForm
161      !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
162         DO_2D( 0, 0, 0, 0 )
163            pr3u(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
164               &                    + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
165            pr3v(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  &
166               &                    + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
167         END_2D
168!!st      ELSE                                         !- Flux Form   (simple averaging)
169#else
170         DO_2D( 0, 0, 0, 0 )
171            pr3u(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji+1,jj  )  ) * r1_hu_0(ji,jj)
172            pr3v(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji  ,jj+1)  ) * r1_hv_0(ji,jj)
173         END_2D
174!!st      ENDIF
175#endif         
176      !
177      IF( .NOT.PRESENT( pr3f ) ) THEN              !- lbc on ratio at u-, v-points only
178         CALL lbc_lnk_multi( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp )
179         !
180         !
181      ELSE                                   !==  ratio at f-point  ==!
182         !
183!!st         IF( ln_dynadv_vec )   THEN                !- Vector Form   (thickness weighted averaging)
184#if ! defined key_qcoTest_FluxForm
185         !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
186
187            DO_2D( 0, 0, 0, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
188               pr3f(ji,jj) = 0.25_wp * (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  &
189                  &                     + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  &
190                  &                     + e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)  &
191                  &                     + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
192            END_2D
193!!st         ELSE                                      !- Flux Form   (simple averaging)
194#else
195            DO_2D( 0, 0, 0, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
196               pr3f(ji,jj) = 0.25_wp * (  pssh(ji,jj  ) + pssh(ji+1,jj  )  &
197                  &                     + pssh(ji,jj+1) + pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj)
198            END_2D
199!!st         ENDIF
200#endif
201         !                                                 ! lbc on ratio at u-,v-,f-points
202         CALL lbc_lnk_multi( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp )
203         !
204      ENDIF
205      !
206   END SUBROUTINE dom_qco_r3c
207
208
209   SUBROUTINE qco_rst_read( kt, Kbb, Kmm )
210      !!---------------------------------------------------------------------
211      !!                   ***  ROUTINE qco_rst_read  ***
212      !!
213      !! ** Purpose :   Read ssh in restart file
214      !!
215      !! ** Method  :   use of IOM library
216      !!                if the restart does not contain ssh,
217      !!                it is set to the _0 values.
218      !!----------------------------------------------------------------------
219      INTEGER, INTENT(in) ::   kt         ! ocean time-step
220      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
221      !
222      INTEGER ::   ji, jj, jk
223      INTEGER ::   id1, id2     ! local integers
224      !!----------------------------------------------------------------------
225      !
226      IF( ln_rstart ) THEN                   !* Read the restart file
227         CALL rst_read_open                  !  open the restart file if necessary
228         !
229         id1 = iom_varid( numror, 'sshb', ldstop = .FALSE. )
230         id2 = iom_varid( numror, 'sshn', ldstop = .FALSE. )
231         !
232         !                             ! --------- !
233         !                             ! all cases !
234         !                             ! --------- !
235         !
236         IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
237            CALL iom_get( numror, jpdom_auto, 'sshb'   , ssh(:,:,Kbb), ldxios = lrxios    )
238            CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
239            ! needed to restart if land processor not computed
240            IF(lwp) write(numout,*) 'qco_rst_read : ssh(:,:,Kbb) and ssh(:,:,Kmm) found in restart files'
241            !!WHERE ( ssmask(:,:) == 0.0_wp )   !!gm/st ==> sm should not be necessary on ssh while it was required on e3
242            !!   ssh(:,:,Kmm) = 0._wp
243            !!   ssh(:,:,Kbb) = 0._wp
244            !!END WHERE
245            IF( l_1st_euler ) THEN
246               ssh(:,:,Kbb) = ssh(:,:,Kmm)
247            ENDIF
248         ELSE IF( id1 > 0 ) THEN
249            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kmm) not found in restart files'
250            IF(lwp) write(numout,*) 'sshn set equal to sshb.'
251            IF(lwp) write(numout,*) 'neuler is forced to 0'
252            CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb), ldxios = lrxios )
253            ssh(:,:,Kmm) = ssh(:,:,Kbb)
254            l_1st_euler = .TRUE.
255         ELSE IF( id2 > 0 ) THEN
256            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kbb) not found in restart files'
257            IF(lwp) write(numout,*) 'sshb set equal to sshn.'
258            IF(lwp) write(numout,*) 'neuler is forced to 0'
259            CALL iom_get( numror, jpdom_auto, 'sshn', ssh(:,:,Kmm), ldxios = lrxios )
260            ssh(:,:,Kbb) = ssh(:,:,Kmm)
261            l_1st_euler = .TRUE.
262         ELSE
263            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kmm) not found in restart file'
264            IF(lwp) write(numout,*) 'ssh_b and ssh_n set to zero'
265            IF(lwp) write(numout,*) 'neuler is forced to 0'
266            ssh(:,:,:) = 0._wp
267            l_1st_euler = .TRUE.
268         ENDIF
269         !
270      ELSE                                   !* Initialize at "rest"
271         !
272         IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
273            !
274            IF( cn_cfg == 'wad' ) THEN            ! Wetting and drying test case
275               CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
276               ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
277               ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb)
278               uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb)
279               vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb)
280            ELSE                                  ! if not test case
281               ssh(:,:,Kmm) = -ssh_ref
282               ssh(:,:,Kbb) = -ssh_ref
283               !
284               DO_2D( 1, 1, 1, 1 )
285                  IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
286                     ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
287                     ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
288                  ENDIF
289               END_2D
290            ENDIF
291            !
292            DO ji = 1, jpi
293               DO jj = 1, jpj
294                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
295                    CALL ctl_stop( 'qco_rst_read: ht_0 must be positive at potentially wet points' )
296                  ENDIF
297               END DO
298            END DO
299            !
300         ELSE
301            !
302            ! Just to read set ssh in fact, called latter once vertical grid
303            ! is set up:
304!           CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
305!           !
306            ssh(:,:,:) = 0._wp
307            !
308         ENDIF           ! end of ll_wd edits
309         !
310      ENDIF
311      !
312   END SUBROUTINE qco_rst_read
313
314   
315   SUBROUTINE qco_rst_read2( kt, Kbb, Kmm )
316      !!---------------------------------------------------------------------
317      !!                   ***  ROUTINE qco_rst_read  ***
318      !!
319      !! ** Purpose :   Read ssh in restart file
320      !!
321      !! ** Method  :   use of IOM library
322      !!                if the restart does not contain ssh,
323      !!                it is set to the _0 values.
324      !!----------------------------------------------------------------------
325      INTEGER, INTENT(in) ::   kt         ! ocean time-step
326      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
327      !
328      INTEGER ::   ji, jj, jk
329      INTEGER ::   id1, id2     ! local integers
330      !!----------------------------------------------------------------------
331      !
332      IF( ln_rstart ) THEN                   !* Read the restart file
333         CALL rst_read_open                  !  open the restart file if necessary
334         !
335         id1 = iom_varid( numror, 'sshb', ldstop = .FALSE. )
336         id2 = iom_varid( numror, 'sshn', ldstop = .FALSE. )
337         !
338         !                             ! --------- !
339         !                             ! all cases !
340         !                             ! --------- !
341         !
342         IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
343            CALL iom_get( numror, jpdom_auto, 'sshb'   , ssh(:,:,Kbb), ldxios = lrxios    )
344            CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
345            CALL iom_get( numror, jpdom_auto, 'r3tb'   , r3t(:,:,Kbb), ldxios = lrxios    )
346            CALL iom_get( numror, jpdom_auto, 'r3tn'   , r3t(:,:,Kmm), ldxios = lrxios    )
347            CALL iom_get( numror, jpdom_auto, 'r3ub'   , r3u(:,:,Kbb), ldxios = lrxios, cd_type = 'U'    )
348            CALL iom_get( numror, jpdom_auto, 'r3un'   , r3u(:,:,Kmm), ldxios = lrxios, cd_type = 'U'    )
349            CALL iom_get( numror, jpdom_auto, 'r3vb'   , r3v(:,:,Kbb), ldxios = lrxios, cd_type = 'V'    )
350            CALL iom_get( numror, jpdom_auto, 'r3vn'   , r3v(:,:,Kmm), ldxios = lrxios, cd_type = 'V'    )
351            CALL iom_get( numror, jpdom_auto, 'r3f'    , r3f(:,:)    , ldxios = lrxios, cd_type = 'F'    )
352           
353            ! needed to restart if land processor not computed
354            IF(lwp) write(numout,*) 'qco_rst_read : ssh(:,:,Kbb) and ssh(:,:,Kmm) found in restart files'
355            !!WHERE ( ssmask(:,:) == 0.0_wp )   !!gm/st ==> sm should not be necessary on ssh while it was required on e3
356            !!   ssh(:,:,Kmm) = 0._wp
357            !!   ssh(:,:,Kbb) = 0._wp
358            !!END WHERE
359            IF( l_1st_euler ) THEN
360               ssh(:,:,Kbb) = ssh(:,:,Kmm)
361            ENDIF
362         ELSE IF( id1 > 0 ) THEN
363            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kmm) not found in restart files'
364            IF(lwp) write(numout,*) 'sshn set equal to sshb.'
365            IF(lwp) write(numout,*) 'neuler is forced to 0'
366            CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb), ldxios = lrxios )
367            ssh(:,:,Kmm) = ssh(:,:,Kbb)
368            l_1st_euler = .TRUE.
369         ELSE IF( id2 > 0 ) THEN
370            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kbb) not found in restart files'
371            IF(lwp) write(numout,*) 'sshb set equal to sshn.'
372            IF(lwp) write(numout,*) 'neuler is forced to 0'
373            CALL iom_get( numror, jpdom_auto, 'sshn', ssh(:,:,Kmm), ldxios = lrxios )
374            ssh(:,:,Kbb) = ssh(:,:,Kmm)
375            l_1st_euler = .TRUE.
376         ELSE
377            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kmm) not found in restart file'
378            IF(lwp) write(numout,*) 'ssh_b and ssh_n set to zero'
379            IF(lwp) write(numout,*) 'neuler is forced to 0'
380            ssh(:,:,:) = 0._wp
381            l_1st_euler = .TRUE.
382         ENDIF
383         !
384      ELSE                                   !* Initialize at "rest"
385         !
386         IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
387            !
388            IF( cn_cfg == 'wad' ) THEN            ! Wetting and drying test case
389               CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
390               ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
391               ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb)
392               uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb)
393               vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb)
394            ELSE                                  ! if not test case
395               ssh(:,:,Kmm) = -ssh_ref
396               ssh(:,:,Kbb) = -ssh_ref
397               !
398               DO_2D( 1, 1, 1, 1 )
399                  IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
400                     ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
401                     ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
402                  ENDIF
403               END_2D
404            ENDIF
405            !
406            DO ji = 1, jpi
407               DO jj = 1, jpj
408                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
409                    CALL ctl_stop( 'qco_rst_read: ht_0 must be positive at potentially wet points' )
410                  ENDIF
411               END DO
412            END DO
413            !
414         ELSE
415            !
416            ! Just to read set ssh in fact, called latter once vertical grid
417            ! is set up:
418            CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
419            !
420            ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb) 
421!            ssh(:,:,:) = 0._wp
422!            r3t(:,:,:) = 0._wp
423!            r3u(:,:,:) = 0._wp
424!            r3v(:,:,:) = 0._wp
425!            r3f(:,:  ) = 0._wp
426            !
427         ENDIF           ! end of ll_wd edits
428         !
429      ENDIF
430      !
431   END SUBROUTINE qco_rst_read2
432
433
434   SUBROUTINE qco_ctl
435      !!---------------------------------------------------------------------
436      !!                  ***  ROUTINE qco_ctl  ***
437      !!
438      !! ** Purpose :   Control the consistency between namelist options
439      !!                for vertical coordinate
440      !!----------------------------------------------------------------------
441      INTEGER ::   ioptio, ios
442      !!
443      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
444         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
445         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
446      !!----------------------------------------------------------------------
447      !
448      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
449901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
450      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
451902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
452      IF(lwm) WRITE ( numond, nam_vvl )
453      !
454      IF(lwp) THEN                    ! Namelist print
455         WRITE(numout,*)
456         WRITE(numout,*) 'qco_ctl : choice/control of the variable vertical coordinate'
457         WRITE(numout,*) '~~~~~~~~'
458         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
459         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
460         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
461         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
462         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
463         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
464         WRITE(numout,*) '      !'
465         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
466         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
467         IF( ln_vvl_ztilde_as_zstar ) THEN
468            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
469            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
470            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
471            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
472            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
473            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rn_Dt'
474         ELSE
475            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
476            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
477         ENDIF
478         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
479      ENDIF
480      !
481      ioptio = 0                      ! Parameter control
482      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
483      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
484      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
485      IF( ln_vvl_layer           )   ioptio = ioptio + 1
486      !
487      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
488      !
489      IF(lwp) THEN                   ! Print the choice
490         WRITE(numout,*)
491         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
492         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
493         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
494         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
495      ENDIF
496      !
497#if defined key_agrif
498      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
499#endif
500      !
501   END SUBROUTINE qco_ctl
502
503   !!======================================================================
504END MODULE domqco
Note: See TracBrowser for help on using the repository browser.