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_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/DYN/dynspg.F90 @ 11671

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

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.
  • 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 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 )
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      !
75      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
76      REAL(wp) ::   z2dt, zg_2, zintp, zgrau0r, zld   ! local scalars
77      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zpice
78      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
79      !!----------------------------------------------------------------------
80      !
81      IF( ln_timing )   CALL timing_start('dyn_spg')
82      !
83      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
84         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
85         ztrdu(:,:,:) = ua(:,:,:)
86         ztrdv(:,:,:) = va(:,:,:)
87      ENDIF
88      !
89      IF(      ln_apr_dyn                                                &   ! atmos. pressure
90         .OR.  ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. ln_tide) )   &   ! tide potential (no time slitting)
91         .OR.  ln_ice_embd ) THEN                                            ! embedded sea-ice
92         !
93         DO jj = 2, jpjm1
94            DO ji = fs_2, fs_jpim1   ! vector opt.
95               spgu(ji,jj) = 0._wp
96               spgv(ji,jj) = 0._wp
97            END DO
98         END DO         
99         !
100         IF( ln_apr_dyn .AND. .NOT.ln_dynspg_ts ) THEN   !==  Atmospheric pressure gradient (added later in time-split case) ==!
101            zg_2 = grav * 0.5
102            DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh
103               DO ji = fs_2, fs_jpim1   ! vector opt.
104                  spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
105                     &                                + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
106                  spgv(ji,jj) = spgv(ji,jj) + zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    &
107                     &                                + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
108               END DO
109            END DO
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            CALL upd_tide( kt )                      ! update tide potential
116            !
117            DO jj = 2, jpjm1                         ! add tide potential forcing
118               DO ji = fs_2, fs_jpim1   ! vector opt.
119                  spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
120                  spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
121               END DO
122            END DO
123            !
124            IF (ln_scal_load) THEN
125               zld = rn_scal_load * grav
126               DO jj = 2, jpjm1                    ! add scalar approximation for load potential
127                  DO ji = fs_2, fs_jpim1   ! vector opt.
128                     spgu(ji,jj) = spgu(ji,jj) + zld * ( sshn(ji+1,jj) - sshn(ji,jj) ) * r1_e1u(ji,jj)
129                     spgv(ji,jj) = spgv(ji,jj) + zld * ( sshn(ji,jj+1) - sshn(ji,jj) ) * r1_e2v(ji,jj)
130                  END DO
131               END DO
132            ENDIF
133         ENDIF
134         !
135         IF( ln_ice_embd ) THEN              !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
136            ALLOCATE( zpice(jpi,jpj) )
137            zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
138            zgrau0r     = - grav * r1_rau0
139            zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r
140            DO jj = 2, jpjm1
141               DO ji = fs_2, fs_jpim1   ! vector opt.
142                  spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
143                  spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
144               END DO
145            END DO
146            DEALLOCATE( zpice )         
147         ENDIF
148         !
149         DO jk = 1, jpkm1                    !== Add all terms to the general trend
150            DO jj = 2, jpjm1
151               DO ji = fs_2, fs_jpim1   ! vector opt.
152                  ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
153                  va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
154               END DO
155            END DO
156         END DO   
157         !
158!!gm add here a call to dyn_trd for ice pressure gradient, the surf pressure trends ????
159         !   
160      ENDIF
161      !
162      SELECT CASE ( nspg )                   !== surface pressure gradient computed and add to the general trend ==!
163      CASE ( np_EXP )   ;   CALL dyn_spg_exp( kt )              ! explicit
164      CASE ( np_TS  )   ;   CALL dyn_spg_ts ( kt )              ! time-splitting
165      END SELECT
166      !                   
167      IF( l_trddyn )   THEN                  ! save the surface pressure gradient trends for further diagnostics
168         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
169         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
170         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt )
171         DEALLOCATE( ztrdu , ztrdv ) 
172      ENDIF
173      !                                      ! print mean trends (used for debugging)
174      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' spg  - Ua: ', mask1=umask, &
175         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
176      !
177      IF( ln_timing )   CALL timing_stop('dyn_spg')
178      !
179   END SUBROUTINE dyn_spg
180
181
182   SUBROUTINE dyn_spg_init
183      !!---------------------------------------------------------------------
184      !!                  ***  ROUTINE dyn_spg_init  ***
185      !!               
186      !! ** Purpose :   Control the consistency between namelist options for
187      !!              surface pressure gradient schemes
188      !!----------------------------------------------------------------------
189      INTEGER ::   ioptio, ios   ! local integers
190      !
191      NAMELIST/namdyn_spg/ ln_dynspg_exp       , ln_dynspg_ts,   &
192      &                    ln_bt_fw, ln_bt_av  , ln_bt_auto  ,   &
193      &                    nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha
194      !!----------------------------------------------------------------------
195      !
196      IF(lwp) THEN
197         WRITE(numout,*)
198         WRITE(numout,*) 'dyn_spg_init : choice of the surface pressure gradient scheme'
199         WRITE(numout,*) '~~~~~~~~~~~~'
200      ENDIF
201      !
202      READ  ( numnam_ref, namdyn_spg, IOSTAT = ios, ERR = 901)
203901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_spg in reference namelist' )
204      !
205      READ  ( numnam_cfg, namdyn_spg, IOSTAT = ios, ERR = 902 )
206902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist' )
207      IF(lwm) WRITE ( numond, namdyn_spg )
208      !
209      IF(lwp) THEN             ! Namelist print
210         WRITE(numout,*) '   Namelist : namdyn_spg                    '
211         WRITE(numout,*) '      Explicit free surface                  ln_dynspg_exp = ', ln_dynspg_exp
212         WRITE(numout,*) '      Free surface with time splitting       ln_dynspg_ts  = ', ln_dynspg_ts
213      ENDIF
214      !                          ! Control of surface pressure gradient scheme options
215                                     nspg =  np_NO    ;   ioptio = 0
216      IF( ln_dynspg_exp ) THEN   ;   nspg =  np_EXP   ;   ioptio = ioptio + 1   ;   ENDIF
217      IF( ln_dynspg_ts  ) THEN   ;   nspg =  np_TS    ;   ioptio = ioptio + 1   ;   ENDIF
218      !
219      IF( ioptio  > 1 )   CALL ctl_stop( 'Choose only one surface pressure gradient scheme' )
220      IF( ioptio == 0 )   CALL ctl_warn( 'NO surface pressure gradient trend in momentum Eqs.' )
221      IF( ln_dynspg_exp .AND. ln_isfcav )   &
222           &   CALL ctl_stop( ' dynspg_exp not tested with ice shelf cavity ' )
223      !
224      IF(lwp) THEN
225         WRITE(numout,*)
226         IF( nspg == np_EXP )   WRITE(numout,*) '   ==>>>   explicit free surface'
227         IF( nspg == np_TS  )   WRITE(numout,*) '   ==>>>   free surface with time splitting scheme'
228         IF( nspg == np_NO  )   WRITE(numout,*) '   ==>>>   No surface surface pressure gradient trend in momentum Eqs.'
229      ENDIF
230      !
231      IF( nspg == np_TS ) THEN   ! split-explicit scheme initialisation
232         CALL dyn_spg_ts_init          ! do it first: set nn_baro used to allocate some arrays later on
233      ENDIF
234      !
235   END SUBROUTINE dyn_spg_init
236
237  !!======================================================================
238END MODULE dynspg
Note: See TracBrowser for help on using the repository browser.