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.
dynspg.F90 in branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
File size: 12.8 KB
Line 
1MODULE dynspg
2   !!======================================================================
3   !!                       ***  MODULE  dynspg  ***
4   !! Ocean dynamics:  surface pressure gradient control
5   !!======================================================================
6   !! History :  1.0  ! 2005-12  (C. Talandier, G. Madec, V. Garnier)  Original code
7   !!            3.2  ! 2009-07  (R. Benshila)  Suppression of rigid-lid option
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   dyn_spg     : update the dynamics trend with the lateral diffusion
12   !!   dyn_spg_ctl : initialization, namelist read, and parameters control
13   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and tracers variables
15   USE dom_oce        ! ocean space and time domain variables
16   USE c1d            ! 1D vertical configuration
17   USE phycst         ! physical constants
18   USE sbc_oce        ! surface boundary condition: ocean
19   USE sbcapr         ! surface boundary condition: atmospheric pressure
20   USE dynspg_oce     ! surface pressure gradient variables
21   USE dynspg_exp     ! surface pressure gradient     (dyn_spg_exp routine)
22   USE dynspg_ts      ! surface pressure gradient     (dyn_spg_ts  routine)
23   USE dynspg_flt     ! surface pressure gradient     (dyn_spg_flt routine)
24   USE dynadv         ! dynamics: vector invariant versus flux form
25   USE dynhpg, ONLY: ln_dynhpg_imp
26   USE sbctide
27   USE updtide
28   USE trd_oce        ! trends: ocean variables
29   USE trddyn         ! trend manager: dynamics
30   !
31   USE prtctl         ! Print control                     (prt_ctl routine)
32   USE in_out_manager ! I/O manager
33   USE lib_mpp        ! MPP library
34   USE solver         ! solver initialization
35   USE wrk_nemo       ! Memory Allocation
36   USE timing         ! Timing
37
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   dyn_spg        ! routine called by step module
43   PUBLIC   dyn_spg_init   ! routine called by opa module
44
45   INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_...
46
47   !! * Substitutions
48#  include "domzgr_substitute.h90"
49#  include "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
52   !! $Id$
53   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE dyn_spg( kt, kindic )
58      !!----------------------------------------------------------------------
59      !!                  ***  ROUTINE dyn_spg  ***
60      !!
61      !! ** Purpose :   achieve the momentum time stepping by computing the
62      !!              last trend, the surface pressure gradient including the
63      !!              atmospheric pressure forcing (ln_apr_dyn=T), and performing
64      !!              the Leap-Frog integration.
65      !!gm              In the current version only the filtered solution provide
66      !!gm            the after velocity, in the 2 other (ua,va) are still the trends
67      !!
68      !! ** Method  :   Three schemes:
69      !!              - explicit computation      : the spg is evaluated at now
70      !!              - filtered computation      : the Roulet & madec (2000) technique is used
71      !!              - split-explicit computation: a time splitting technique is used
72      !!
73      !!              ln_apr_dyn=T : the atmospheric pressure forcing is applied
74      !!             as the gradient of the inverse barometer ssh:
75      !!                apgu = - 1/rau0 di[apr] = 0.5*grav di[ssh_ib+ssh_ibb]
76      !!                apgv = - 1/rau0 dj[apr] = 0.5*grav dj[ssh_ib+ssh_ibb]
77      !!             Note that as all external forcing a time averaging over a two rdt
78      !!             period is used to prevent the divergence of odd and even time step.
79      !!----------------------------------------------------------------------
80      !
81      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
82      INTEGER, INTENT(  out) ::   kindic   ! solver flag
83      !
84      INTEGER  ::   ji, jj, jk                             ! dummy loop indices
85      REAL(wp) ::   z2dt, zg_2, zintp, zgrau0r             ! temporary scalar
86      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
87      REAL(wp), POINTER, DIMENSION(:,:)   ::  zpice
88      !!----------------------------------------------------------------------
89      !
90      IF( nn_timing == 1 )  CALL timing_start('dyn_spg')
91      !
92
93!!gm NOTA BENE : the dynspg_exp and dynspg_ts should be modified so that
94!!gm             they return the after velocity, not the trends (as in trazdf_imp...)
95!!gm             In this case, change/simplify dynnxt
96
97
98      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
99         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
100         ztrdu(:,:,:) = ua(:,:,:)
101         ztrdv(:,:,:) = va(:,:,:)
102      ENDIF
103
104      IF(      ln_apr_dyn                                                &   ! atmos. pressure
105         .OR.  ( .NOT.lk_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) )   &   ! tide potential (no time slitting)
106         .OR.  nn_ice_embd == 2  ) THEN                                      ! embedded sea-ice
107         !
108         DO jj = 2, jpjm1
109            DO ji = fs_2, fs_jpim1   ! vector opt.
110               spgu(ji,jj) = 0._wp
111               spgv(ji,jj) = 0._wp
112            END DO
113         END DO         
114         !
115         IF( ln_apr_dyn .AND. (.NOT. lk_dynspg_ts) ) THEN                    !==  Atmospheric pressure gradient (added later in time-split case) ==!
116            zg_2 = grav * 0.5
117            DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh
118               DO ji = fs_2, fs_jpim1   ! vector opt.
119                  spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
120                     &                      + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
121                  spgv(ji,jj) = spgv(ji,jj) + zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    &
122                     &                      + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
123               END DO
124            END DO
125         ENDIF
126         !
127         !                                    !==  tide potential forcing term  ==!
128         IF( .NOT.lk_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case
129            !
130            CALL upd_tide( kt )                      ! update tide potential
131            !
132            DO jj = 2, jpjm1                         ! add tide potential forcing
133               DO ji = fs_2, fs_jpim1   ! vector opt.
134                  spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
135                  spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
136               END DO
137            END DO
138         ENDIF
139         !
140         IF( nn_ice_embd == 2 ) THEN          !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
141            CALL wrk_alloc( jpi, jpj, zpice )
142            !                                           
143            zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
144            zgrau0r     = - grav * r1_rau0
145            zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r
146            DO jj = 2, jpjm1
147               DO ji = fs_2, fs_jpim1   ! vector opt.
148                  spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
149                  spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
150               END DO
151            END DO
152            !
153            CALL wrk_dealloc( jpi, jpj, zpice )         
154         ENDIF
155         !
156         DO jk = 1, jpkm1                     !== Add all terms to the general trend
157            DO jj = 2, jpjm1
158               DO ji = fs_2, fs_jpim1   ! vector opt.
159                  ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
160                  va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
161               END DO
162            END DO
163         END DO   
164         
165!!gm add here a call to dyn_trd for ice pressure gradient, the surf pressure trends ????
166             
167      ENDIF
168
169      SELECT CASE ( nspg )                       ! compute surf. pressure gradient trend and add it to the general trend
170      !                                                     
171      CASE (  0 )   ;   CALL dyn_spg_exp( kt )              ! explicit
172      CASE (  1 )   ;   CALL dyn_spg_ts ( kt )              ! time-splitting
173      CASE (  2 )   ;   CALL dyn_spg_flt( kt, kindic )      ! filtered
174      !                                                   
175      END SELECT
176      !                   
177      IF( l_trddyn )   THEN                      ! save the surface pressure gradient trends for further diagnostics
178         SELECT CASE ( nspg )
179         CASE ( 0, 1 )
180            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
181            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
182         CASE( 2 )
183            z2dt = 2. * rdt
184            IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
185            ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / z2dt - ztrdu(:,:,:)
186            ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:)
187         END SELECT
188         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt )
189         !
190         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
191      ENDIF
192      !                                          ! print mean trends (used for debugging)
193      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' spg  - Ua: ', mask1=umask, &
194         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
195      !
196      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg')
197      !
198   END SUBROUTINE dyn_spg
199
200
201   SUBROUTINE dyn_spg_init
202      !!---------------------------------------------------------------------
203      !!                  ***  ROUTINE dyn_spg_init  ***
204      !!               
205      !! ** Purpose :   Control the consistency between cpp options for
206      !!              surface pressure gradient schemes
207      !!----------------------------------------------------------------------
208      INTEGER ::   ioptio
209      !!----------------------------------------------------------------------
210      !
211      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_init')
212      !
213      IF(lwp) THEN             ! Control print
214         WRITE(numout,*)
215         WRITE(numout,*) 'dyn_spg_init : choice of the surface pressure gradient scheme'
216         WRITE(numout,*) '~~~~~~~~~~~'
217         WRITE(numout,*) '     Explicit free surface                  lk_dynspg_exp = ', lk_dynspg_exp
218         WRITE(numout,*) '     Free surface with time splitting       lk_dynspg_ts  = ', lk_dynspg_ts
219         WRITE(numout,*) '     Filtered free surface cst volume       lk_dynspg_flt = ', lk_dynspg_flt
220      ENDIF
221
222      IF( lk_dynspg_ts ) CALL dyn_spg_ts_init( nit000 )
223      ! (do it now, to set nn_baro, used to allocate some arrays later on)
224      !                        ! allocate dyn_spg arrays
225      IF( lk_dynspg_ts ) THEN
226         IF( dynspg_oce_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_oce arrays')
227         IF( dyn_spg_ts_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays')
228         IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' )
229      ENDIF
230
231      !                        ! Control of surface pressure gradient scheme options
232      ioptio = 0
233      IF(lk_dynspg_exp)   ioptio = ioptio + 1
234      IF(lk_dynspg_ts )   ioptio = ioptio + 1
235      IF(lk_dynspg_flt)   ioptio = ioptio + 1
236      !
237      IF(  ioptio > 1  .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   &
238           &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' )
239      IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav )   &
240           &   CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' )
241      !
242      IF( lk_dynspg_exp)   nspg =  0
243      IF( lk_dynspg_ts )   nspg =  1
244      IF( lk_dynspg_flt)   nspg =  2
245      !
246      IF(lwp) THEN
247         WRITE(numout,*)
248         IF( nspg ==  0 )   WRITE(numout,*) '     explicit free surface'
249         IF( nspg ==  1 )   WRITE(numout,*) '     free surface with time splitting scheme'
250         IF( nspg ==  2 )   WRITE(numout,*) '     filtered free surface'
251      ENDIF
252
253#if defined key_dynspg_flt
254      CALL solver_init( nit000 )   ! Elliptic solver initialisation
255#endif
256      !               ! Control of hydrostatic pressure choice
257      IF( lk_dynspg_ts .AND. ln_dynhpg_imp )   CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' )
258      !
259      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_init')
260      !
261   END SUBROUTINE dyn_spg_init
262
263  !!======================================================================
264END MODULE dynspg
Note: See TracBrowser for help on using the repository browser.