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/trunk/src/OCE/DYN – NEMO

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

Last change on this file since 12377 was 12377, checked in by acc, 4 years 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
RevLine 
[358]1MODULE dynspg
2   !!======================================================================
3   !!                       ***  MODULE  dynspg  ***
4   !! Ocean dynamics:  surface pressure gradient control
5   !!======================================================================
[1566]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
[503]8   !!----------------------------------------------------------------------
[358]9
10   !!----------------------------------------------------------------------
[5930]11   !!   dyn_spg     : update the dynamics trend with surface pressure gradient
12   !!   dyn_spg_init: initialization, namelist read, and parameters control
[358]13   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and tracers variables
15   USE dom_oce        ! ocean space and time domain variables
[4245]16   USE c1d            ! 1D vertical configuration
[2715]17   USE phycst         ! physical constants
[2528]18   USE sbc_oce        ! surface boundary condition: ocean
[9019]19   USE sbc_ice , ONLY : snwice_mass, snwice_mass_b
[2528]20   USE sbcapr         ! surface boundary condition: atmospheric pressure
[358]21   USE dynspg_exp     ! surface pressure gradient     (dyn_spg_exp routine)
22   USE dynspg_ts      ! surface pressure gradient     (dyn_spg_ts  routine)
[12377]23   USE tide_mod       !
[4990]24   USE trd_oce        ! trends: ocean variables
25   USE trddyn         ! trend manager: dynamics
26   !
[358]27   USE prtctl         ! Print control                     (prt_ctl routine)
28   USE in_out_manager ! I/O manager
[2715]29   USE lib_mpp        ! MPP library
[4990]30   USE timing         ! Timing
[358]31
32   IMPLICIT NONE
33   PRIVATE
34
[2528]35   PUBLIC   dyn_spg        ! routine called by step module
36   PUBLIC   dyn_spg_init   ! routine called by opa module
[358]37
[6140]38   INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_...
[358]39
[6140]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
[12377]44   !
45   REAL(wp) ::   zt0step !   Time of day at the beginning of the time step
[6140]46
[358]47   !! * Substitutions
[12377]48#  include "do_loop_substitute.h90"
[358]49   !!----------------------------------------------------------------------
[9598]50   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]51   !! $Id$
[10068]52   !! Software governed by the CeCILL license (see ./LICENSE)
[358]53   !!----------------------------------------------------------------------
54CONTAINS
55
[12377]56   SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa )
[358]57      !!----------------------------------------------------------------------
58      !!                  ***  ROUTINE dyn_spg  ***
59      !!
[6140]60      !! ** Purpose :   compute surface pressure gradient including the
61      !!              atmospheric pressure forcing (ln_apr_dyn=T).
[1566]62      !!
[5930]63      !! ** Method  :   Two schemes:
[6140]64      !!              - explicit       : the spg is evaluated at now
65      !!              - split-explicit : a time splitting technique is used
[1566]66      !!
[2528]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.
[358]73      !!----------------------------------------------------------------------
[12377]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
[2715]78      !
[9019]79      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
[9023]80      REAL(wp) ::   z2dt, zg_2, zintp, zgrau0r, zld   ! local scalars
[9019]81      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zpice
82      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
[358]83      !!----------------------------------------------------------------------
[3294]84      !
[9019]85      IF( ln_timing )   CALL timing_start('dyn_spg')
[3294]86      !
[358]87      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
[9019]88         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
[12377]89         ztrdu(:,:,:) = puu(:,:,:,Krhs)
90         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
[358]91      ENDIF
[6140]92      !
[4292]93      IF(      ln_apr_dyn                                                &   ! atmos. pressure
[7646]94         .OR.  ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. ln_tide) )   &   ! tide potential (no time slitting)
[9019]95         .OR.  ln_ice_embd ) THEN                                            ! embedded sea-ice
[4292]96         !
[12377]97         DO_2D_00_00
98            spgu(ji,jj) = 0._wp
99            spgv(ji,jj) = 0._wp
100         END_2D
[4292]101         !
[6140]102         IF( ln_apr_dyn .AND. .NOT.ln_dynspg_ts ) THEN   !==  Atmospheric pressure gradient (added later in time-split case) ==!
[4292]103            zg_2 = grav * 0.5
[12377]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
[4292]110         ENDIF
111         !
112         !                                    !==  tide potential forcing term  ==!
[7646]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
[4292]114            !
[12377]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)
[4292]118            !
[12377]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
[9023]123            !
124            IF (ln_scal_load) THEN
125               zld = rn_scal_load * grav
[12377]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
[9023]130            ENDIF
[4292]131         ENDIF
132         !
[9019]133         IF( ln_ice_embd ) THEN              !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
134            ALLOCATE( zpice(jpi,jpj) )
[4292]135            zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
136            zgrau0r     = - grav * r1_rau0
[7753]137            zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r
[12377]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
[9019]142            DEALLOCATE( zpice )         
[4292]143         ENDIF
144         !
[12377]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
[6140]149         !
[4990]150!!gm add here a call to dyn_trd for ice pressure gradient, the surf pressure trends ????
[6140]151         !   
[3625]152      ENDIF
[6140]153      !
154      SELECT CASE ( nspg )                   !== surface pressure gradient computed and add to the general trend ==!
[12377]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
[358]157      END SELECT
[503]158      !                   
[6140]159      IF( l_trddyn )   THEN                  ! save the surface pressure gradient trends for further diagnostics
[12377]160         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
161         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
162         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt, Kmm )
[9019]163         DEALLOCATE( ztrdu , ztrdv ) 
[358]164      ENDIF
[6140]165      !                                      ! print mean trends (used for debugging)
[12377]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' )
[503]168      !
[9019]169      IF( ln_timing )   CALL timing_stop('dyn_spg')
[2715]170      !
[358]171   END SUBROUTINE dyn_spg
172
173
[2528]174   SUBROUTINE dyn_spg_init
[358]175      !!---------------------------------------------------------------------
[2528]176      !!                  ***  ROUTINE dyn_spg_init  ***
[358]177      !!               
[5930]178      !! ** Purpose :   Control the consistency between namelist options for
[1566]179      !!              surface pressure gradient schemes
[358]180      !!----------------------------------------------------------------------
[6140]181      INTEGER ::   ioptio, ios   ! local integers
[5930]182      !
[6140]183      NAMELIST/namdyn_spg/ ln_dynspg_exp       , ln_dynspg_ts,   &
184      &                    ln_bt_fw, ln_bt_av  , ln_bt_auto  ,   &
[9023]185      &                    nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha
[358]186      !!----------------------------------------------------------------------
[3294]187      !
[9169]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      !
[5930]194      READ  ( numnam_ref, namdyn_spg, IOSTAT = ios, ERR = 901)
[11536]195901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_spg in reference namelist' )
[6140]196      !
[5930]197      READ  ( numnam_cfg, namdyn_spg, IOSTAT = ios, ERR = 902 )
[11536]198902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist' )
[5930]199      IF(lwm) WRITE ( numond, namdyn_spg )
200      !
201      IF(lwp) THEN             ! Namelist print
[9169]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
[358]205      ENDIF
[6140]206      !                          ! Control of surface pressure gradient scheme options
[6981]207                                     nspg =  np_NO    ;   ioptio = 0
[6140]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
[1566]210      !
[6140]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 ' )
[1566]215      !
216      IF(lwp) THEN
[358]217         WRITE(numout,*)
[9190]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.'
[358]221      ENDIF
[1566]222      !
[6140]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      !
[2528]227   END SUBROUTINE dyn_spg_init
[358]228
229  !!======================================================================
230END MODULE dynspg
Note: See TracBrowser for help on using the repository browser.