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 NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg.F90 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago

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