source: NEMO/trunk/src/OCE/DYN/dynspg.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 9 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

  • Property svn:keywords set to Id
File size: 11.9 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 surface pressure gradient
12   !!   dyn_spg_init: 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 sbc_ice , ONLY : snwice_mass, snwice_mass_b
20   USE sbcapr         ! surface boundary condition: atmospheric pressure
21   USE dynspg_exp     ! surface pressure gradient     (dyn_spg_exp routine)
22   USE dynspg_ts      ! surface pressure gradient     (dyn_spg_ts  routine)
23   USE tide_mod       !
24   USE trd_oce        ! trends: ocean variables
25   USE trddyn         ! trend manager: dynamics
26   !
27   USE prtctl         ! Print control                     (prt_ctl routine)
28   USE in_out_manager ! I/O manager
29   USE lib_mpp        ! MPP library
30   USE timing         ! Timing
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   dyn_spg        ! routine called by step module
36   PUBLIC   dyn_spg_init   ! routine called by opa module
37
38   INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_...
39
40   !                       ! Parameter to control the surface pressure gradient scheme
41   INTEGER, PARAMETER ::   np_TS  = 1   ! split-explicit time stepping (Time-Splitting)
42   INTEGER, PARAMETER ::   np_EXP = 0   !       explicit time stepping
43   INTEGER, PARAMETER ::   np_NO  =-1   ! no surface pressure gradient, no scheme
44   !
45   REAL(wp) ::   zt0step !   Time of day at the beginning of the time step
46
47   !! * Substitutions
48#  include "do_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
51   !! $Id$
52   !! Software governed by the CeCILL license (see ./LICENSE)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa )
57      !!----------------------------------------------------------------------
58      !!                  ***  ROUTINE dyn_spg  ***
59      !!
60      !! ** Purpose :   compute surface pressure gradient including the
61      !!              atmospheric pressure forcing (ln_apr_dyn=T).
62      !!
63      !! ** Method  :   Two schemes:
64      !!              - explicit       : the spg is evaluated at now
65      !!              - split-explicit : a time splitting technique is used
66      !!
67      !!              ln_apr_dyn=T : the atmospheric pressure forcing is applied
68      !!             as the gradient of the inverse barometer ssh:
69      !!                apgu = - 1/rau0 di[apr] = 0.5*grav di[ssh_ib+ssh_ibb]
70      !!                apgv = - 1/rau0 dj[apr] = 0.5*grav dj[ssh_ib+ssh_ibb]
71      !!             Note that as all external forcing a time averaging over a two rdt
72      !!             period is used to prevent the divergence of odd and even time step.
73      !!----------------------------------------------------------------------
74      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index
75      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices
76      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation
77      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(inout) ::  pssh, puu_b, pvv_b  ! SSH and barotropic velocities at main time levels
78      !
79      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
80      REAL(wp) ::   z2dt, zg_2, zintp, zgrau0r, zld   ! local scalars
81      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zpice
82      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
83      !!----------------------------------------------------------------------
84      !
85      IF( ln_timing )   CALL timing_start('dyn_spg')
86      !
87      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
88         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
89         ztrdu(:,:,:) = puu(:,:,:,Krhs)
90         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
91      ENDIF
92      !
93      IF(      ln_apr_dyn                                                &   ! atmos. pressure
94         .OR.  ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. ln_tide) )   &   ! tide potential (no time slitting)
95         .OR.  ln_ice_embd ) THEN                                            ! embedded sea-ice
96         !
97         DO_2D_00_00
98            spgu(ji,jj) = 0._wp
99            spgv(ji,jj) = 0._wp
100         END_2D
101         !
102         IF( ln_apr_dyn .AND. .NOT.ln_dynspg_ts ) THEN   !==  Atmospheric pressure gradient (added later in time-split case) ==!
103            zg_2 = grav * 0.5
104            DO_2D_00_00
105               spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
106                  &                                + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
107               spgv(ji,jj) = spgv(ji,jj) + zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    &
108                  &                                + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
109            END_2D
110         ENDIF
111         !
112         !                                    !==  tide potential forcing term  ==!
113         IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. ln_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case
114            !
115            ! Update tide potential at the beginning of current time step
116            zt0step = REAL(nsec_day, wp)-0.5_wp*rdt
117            CALL upd_tide(zt0step, Kmm)
118            !
119            DO_2D_00_00
120               spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
121               spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
122            END_2D
123            !
124            IF (ln_scal_load) THEN
125               zld = rn_scal_load * grav
126               DO_2D_00_00
127                  spgu(ji,jj) = spgu(ji,jj) + zld * ( pssh(ji+1,jj,Kmm) - pssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
128                  spgv(ji,jj) = spgv(ji,jj) + zld * ( pssh(ji,jj+1,Kmm) - pssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
129               END_2D
130            ENDIF
131         ENDIF
132         !
133         IF( ln_ice_embd ) THEN              !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
134            ALLOCATE( zpice(jpi,jpj) )
135            zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
136            zgrau0r     = - grav * r1_rau0
137            zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r
138            DO_2D_00_00
139               spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
140               spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
141            END_2D
142            DEALLOCATE( zpice )         
143         ENDIF
144         !
145         DO_3D_00_00( 1, jpkm1 )
146            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj)
147            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj)
148         END_3D
149         !
150!!gm add here a call to dyn_trd for ice pressure gradient, the surf pressure trends ????
151         !   
152      ENDIF
153      !
154      SELECT CASE ( nspg )                   !== surface pressure gradient computed and add to the general trend ==!
155      CASE ( np_EXP )   ;   CALL dyn_spg_exp( kt,      Kmm,       puu, pvv, Krhs )                    ! explicit
156      CASE ( np_TS  )   ;   CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) ! time-splitting
157      END SELECT
158      !                   
159      IF( l_trddyn )   THEN                  ! save the surface pressure gradient trends for further diagnostics
160         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
161         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
162         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt, Kmm )
163         DEALLOCATE( ztrdu , ztrdv ) 
164      ENDIF
165      !                                      ! print mean trends (used for debugging)
166      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' spg  - Ua: ', mask1=umask, &
167         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
168      !
169      IF( ln_timing )   CALL timing_stop('dyn_spg')
170      !
171   END SUBROUTINE dyn_spg
172
173
174   SUBROUTINE dyn_spg_init
175      !!---------------------------------------------------------------------
176      !!                  ***  ROUTINE dyn_spg_init  ***
177      !!               
178      !! ** Purpose :   Control the consistency between namelist options for
179      !!              surface pressure gradient schemes
180      !!----------------------------------------------------------------------
181      INTEGER ::   ioptio, ios   ! local integers
182      !
183      NAMELIST/namdyn_spg/ ln_dynspg_exp       , ln_dynspg_ts,   &
184      &                    ln_bt_fw, ln_bt_av  , ln_bt_auto  ,   &
185      &                    nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha
186      !!----------------------------------------------------------------------
187      !
188      IF(lwp) THEN
189         WRITE(numout,*)
190         WRITE(numout,*) 'dyn_spg_init : choice of the surface pressure gradient scheme'
191         WRITE(numout,*) '~~~~~~~~~~~~'
192      ENDIF
193      !
194      READ  ( numnam_ref, namdyn_spg, IOSTAT = ios, ERR = 901)
195901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_spg in reference namelist' )
196      !
197      READ  ( numnam_cfg, namdyn_spg, IOSTAT = ios, ERR = 902 )
198902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist' )
199      IF(lwm) WRITE ( numond, namdyn_spg )
200      !
201      IF(lwp) THEN             ! Namelist print
202         WRITE(numout,*) '   Namelist : namdyn_spg                    '
203         WRITE(numout,*) '      Explicit free surface                  ln_dynspg_exp = ', ln_dynspg_exp
204         WRITE(numout,*) '      Free surface with time splitting       ln_dynspg_ts  = ', ln_dynspg_ts
205      ENDIF
206      !                          ! Control of surface pressure gradient scheme options
207                                     nspg =  np_NO    ;   ioptio = 0
208      IF( ln_dynspg_exp ) THEN   ;   nspg =  np_EXP   ;   ioptio = ioptio + 1   ;   ENDIF
209      IF( ln_dynspg_ts  ) THEN   ;   nspg =  np_TS    ;   ioptio = ioptio + 1   ;   ENDIF
210      !
211      IF( ioptio  > 1 )   CALL ctl_stop( 'Choose only one surface pressure gradient scheme' )
212      IF( ioptio == 0 )   CALL ctl_warn( 'NO surface pressure gradient trend in momentum Eqs.' )
213      IF( ln_dynspg_exp .AND. ln_isfcav )   &
214           &   CALL ctl_stop( ' dynspg_exp not tested with ice shelf cavity ' )
215      !
216      IF(lwp) THEN
217         WRITE(numout,*)
218         IF( nspg == np_EXP )   WRITE(numout,*) '   ==>>>   explicit free surface'
219         IF( nspg == np_TS  )   WRITE(numout,*) '   ==>>>   free surface with time splitting scheme'
220         IF( nspg == np_NO  )   WRITE(numout,*) '   ==>>>   No surface surface pressure gradient trend in momentum Eqs.'
221      ENDIF
222      !
223      IF( nspg == np_TS ) THEN   ! split-explicit scheme initialisation
224         CALL dyn_spg_ts_init          ! do it first: set nn_baro used to allocate some arrays later on
225      ENDIF
226      !
227   END SUBROUTINE dyn_spg_init
228
229  !!======================================================================
230END MODULE dynspg
Note: See TracBrowser for help on using the repository browser.