Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

Timestamp:

2010-12-27T18:33:53+01:00 (13 years ago)

Author:

rblod

Message:

Update NEMOGCM from branch nemo_v3_3_beta

File:

• trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90 (modified) (10 diffs, 2 props)

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

@@ -15 +15 @@
    USE dom_oce        ! ocean space and time domain variables
    USE obc_oce        ! ocean open boundary conditions
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE sbcapr         ! surface boundary condition: atmospheric pressure
    USE dynspg_oce     ! surface pressure gradient variables
    USE dynspg_exp     ! surface pressure gradient     (dyn_spg_exp routine)
    USE dynspg_ts      ! surface pressure gradient     (dyn_spg_ts  routine)
    USE dynspg_flt     ! surface pressure gradient     (dyn_spg_flt routine)
+   USE dynadv         ! dynamics: vector invariant versus flux form
    USE trdmod         ! ocean dynamics trends
    USE trdmod_oce     ! ocean variables trends
    USE prtctl         ! Print control                     (prt_ctl routine)
    USE in_out_manager ! I/O manager
+   USE phycst         ! physical constants
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   dyn_spg   ! routine called by step module
+   PUBLIC   dyn_spg        ! routine called by step module
+   PUBLIC   dyn_spg_init   ! routine called by opa module
 
    INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_... 
@@ -37 +42 @@
    !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -47 +52 @@
       !!
       !! ** Purpose :   achieve the momentum time stepping by computing the
-      !!              last trend, the surface pressure gradient, and performing
+      !!              last trend, the surface pressure gradient including the 
+      !!              atmospheric pressure forcing (ln_apr_dyn=T), and performing
       !!              the Leap-Frog integration.
       !!gm              In the current version only the filtered solution provide
@@ -57 +63 @@
       !!              - split-explicit computation: a time splitting technique is used
       !!
+      !!              ln_apr_dyn=T : the atmospheric pressure forcing is applied 
+      !!             as the gradient of the inverse barometer ssh:
+      !!                apgu = - 1/rau0 di[apr] = 0.5*grav di[ssh_ib+ssh_ibb]
+      !!                apgv = - 1/rau0 dj[apr] = 0.5*grav dj[ssh_ib+ssh_ibb]
+      !!             Note that as all external forcing a time averaging over a two rdt
+      !!             period is used to prevent the divergence of odd and even time step.
+      !!
       !! N.B. : When key_esopa is used all the scheme are tested, regardless 
       !!        of the physical meaning of the results. 
@@ -63 +76 @@
       INTEGER, INTENT(  out) ::   kindic   ! solver flag
       !!
-      REAL(wp) ::   z2dt   ! temporary scalar
+      INTEGER  ::   ji, jj, jk                             ! dummy loop indices
+      REAL(wp) ::   z2dt, zg_2                             ! temporary scalar
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrdu, ztrdv   ! 3D workspace
       !!----------------------------------------------------------------------
@@ -74 +88 @@
 
 
-      IF( kt == nit000 )   CALL dyn_spg_ctl      ! initialisation & control of options
-
       IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
          ztrdu(:,:,:) = ua(:,:,:)
          ztrdv(:,:,:) = va(:,:,:)
       ENDIF
+
+      IF( ln_apr_dyn ) THEN                   !==  Atmospheric pressure gradient  ==!
+         zg_2 = grav * 0.5
+         DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               spgu(ji,jj) =  zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
+                  &                   + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) /e1u(ji,jj)
+               spgv(ji,jj) =  zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    &
+                  &                   + ssh_ibb(ji,jj+1) - ssh_ib (ji,jj)  ) /e2v(ji,jj)
+            END DO
+         END DO
+         DO jk = 1, jpkm1                          ! Add the apg to the general trend
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
+                  va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
+               END DO
+            END DO
+         END DO
+      ENDIF
+
 
       SELECT CASE ( nspg )                       ! compute surf. pressure gradient trend and add it to the general trend
@@ -119 +152 @@
 
 
-   SUBROUTINE dyn_spg_ctl
+   SUBROUTINE dyn_spg_init
       !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE dyn_spg_ctl  ***
+      !!                  ***  ROUTINE dyn_spg_init  ***
       !!                
       !! ** Purpose :   Control the consistency between cpp options for 
@@ -131 +164 @@
       IF(lwp) THEN             ! Control print
          WRITE(numout,*)
-         WRITE(numout,*) 'dyn_spg_ctl : choice of the surface pressure gradient scheme'
+         WRITE(numout,*) 'dyn_spg_init : choice of the surface pressure gradient scheme'
          WRITE(numout,*) '~~~~~~~~~~~'
          WRITE(numout,*) '     Explicit free surface                  lk_dynspg_exp = ', lk_dynspg_exp
@@ -164 +197 @@
       !                        ! Control of timestep choice
       IF( lk_dynspg_ts .OR. lk_dynspg_exp ) THEN
-         IF( n_cla == 1 )   &
-           &   CALL ctl_stop( ' Crossland advection not implemented for this free surface formulation ' )
+         IF( nn_cla == 1 )   CALL ctl_stop( 'Crossland advection not implemented for this free surface formulation' )
+      ENDIF
+
+      !                        ! Control of momentum formulation
+      IF( lk_dynspg_ts .AND. lk_vvl ) THEN
+         IF( .NOT.ln_dynadv_vec )   CALL ctl_stop( 'Flux form not implemented for this free surface formulation' )
       ENDIF
 
@@ -178 +215 @@
 #endif
       !
-   END SUBROUTINE dyn_spg_ctl
+   END SUBROUTINE dyn_spg_init
 
   !!======================================================================