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.
Changeset 5989 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90 – NEMO

Ignore:
Timestamp:
2015-12-03T09:10:32+01:00 (8 years ago)
Author:
deazer
Message:

Merging TMB and 25h diagnostics to head of trunk
added brief documentation

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r5260 r5989  
    66   !! History :  1.0  ! 2005-12  (C. Talandier, G. Madec, V. Garnier)  Original code 
    77   !!            3.2  ! 2009-07  (R. Benshila)  Suppression of rigid-lid option 
    8    !!---------------------------------------------------------------------- 
    9  
    10    !!---------------------------------------------------------------------- 
    11    !!   dyn_spg     : update the dynamics trend with the lateral diffusion 
    12    !!   dyn_spg_ctl : initialization, namelist read, and parameters control 
     8   !!            3.7  ! 2015-11  (J. Chanut) Suppression of filtered free surface  
     9   !!---------------------------------------------------------------------- 
     10 
     11   !!---------------------------------------------------------------------- 
     12   !!   dyn_spg     : update the dynamics trend with surface pressure gradient  
     13   !!   dyn_spg_init: initialization, namelist read, and parameters control 
    1314   !!---------------------------------------------------------------------- 
    1415   USE oce            ! ocean dynamics and tracers variables 
     
    1819   USE sbc_oce        ! surface boundary condition: ocean 
    1920   USE sbcapr         ! surface boundary condition: atmospheric pressure 
    20    USE dynspg_oce     ! surface pressure gradient variables 
    2121   USE dynspg_exp     ! surface pressure gradient     (dyn_spg_exp routine) 
    2222   USE dynspg_ts      ! surface pressure gradient     (dyn_spg_ts  routine) 
    23    USE dynspg_flt     ! surface pressure gradient     (dyn_spg_flt routine) 
    24    USE dynadv         ! dynamics: vector invariant versus flux form 
    25    USE dynhpg, ONLY: ln_dynhpg_imp 
    2623   USE sbctide 
    2724   USE updtide 
     
    3229   USE in_out_manager ! I/O manager 
    3330   USE lib_mpp        ! MPP library 
    34    USE solver         ! solver initialization 
    3531   USE wrk_nemo       ! Memory Allocation 
    3632   USE timing         ! Timing 
     
    4339   PUBLIC   dyn_spg_init   ! routine called by opa module 
    4440 
    45    INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_...  
     41   INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from ln_dynspg_...  
    4642 
    4743   !! * Substitutions 
     
    5551CONTAINS 
    5652 
    57    SUBROUTINE dyn_spg( kt, kindic ) 
     53   SUBROUTINE dyn_spg( kt ) 
    5854      !!---------------------------------------------------------------------- 
    5955      !!                  ***  ROUTINE dyn_spg  *** 
     
    6662      !!gm            the after velocity, in the 2 other (ua,va) are still the trends 
    6763      !! 
    68       !! ** Method  :   Three schemes: 
     64      !! ** Method  :   Two schemes: 
    6965      !!              - explicit computation      : the spg is evaluated at now 
    70       !!              - filtered computation      : the Roulet & madec (2000) technique is used 
    7166      !!              - split-explicit computation: a time splitting technique is used 
    7267      !! 
     
    7772      !!             Note that as all external forcing a time averaging over a two rdt 
    7873      !!             period is used to prevent the divergence of odd and even time step. 
    79       !! 
    80       !! N.B. : When key_esopa is used all the scheme are tested, regardless  
    81       !!        of the physical meaning of the results.  
    8274      !!---------------------------------------------------------------------- 
    8375      ! 
    8476      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index 
    85       INTEGER, INTENT(  out) ::   kindic   ! solver flag 
    8677      ! 
    8778      INTEGER  ::   ji, jj, jk                             ! dummy loop indices 
     
    10697 
    10798      IF(      ln_apr_dyn                                                &   ! atmos. pressure 
    108          .OR.  ( .NOT.lk_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) )   &   ! tide potential (no time slitting) 
     99         .OR.  ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) )   &   ! tide potential (no time slitting) 
    109100         .OR.  nn_ice_embd == 2  ) THEN                                      ! embedded sea-ice 
    110101         ! 
     
    116107         END DO          
    117108         ! 
    118          IF( ln_apr_dyn .AND. (.NOT. lk_dynspg_ts) ) THEN                    !==  Atmospheric pressure gradient (added later in time-split case) ==! 
     109         IF( ln_apr_dyn .AND. (.NOT. ln_dynspg_ts) ) THEN                    !==  Atmospheric pressure gradient (added later in time-split case) ==! 
    119110            zg_2 = grav * 0.5 
    120111            DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh 
    121112               DO ji = fs_2, fs_jpim1   ! vector opt. 
    122113                  spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    & 
    123                      &                      + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) /e1u(ji,jj) 
     114                     &                      + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj) 
    124115                  spgv(ji,jj) = spgv(ji,jj) + zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    & 
    125                      &                      + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj) 
     116                     &                      + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj) 
    126117               END DO 
    127118            END DO 
     
    129120         ! 
    130121         !                                    !==  tide potential forcing term  ==! 
    131          IF( .NOT.lk_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case 
     122         IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case 
    132123            ! 
    133124            CALL upd_tide( kt )                      ! update tide potential 
     
    135126            DO jj = 2, jpjm1                         ! add tide potential forcing 
    136127               DO ji = fs_2, fs_jpim1   ! vector opt. 
    137                   spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    138                   spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
     128                  spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 
     129                  spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 
    139130               END DO  
    140131            END DO 
     
    149140            DO jj = 2, jpjm1 
    150141               DO ji = fs_2, fs_jpim1   ! vector opt. 
    151                   spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj) 
    152                   spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj) 
     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) 
    153144               END DO 
    154145            END DO 
     
    174165      CASE (  0 )   ;   CALL dyn_spg_exp( kt )              ! explicit 
    175166      CASE (  1 )   ;   CALL dyn_spg_ts ( kt )              ! time-splitting 
    176       CASE (  2 )   ;   CALL dyn_spg_flt( kt, kindic )      ! filtered 
    177167      !                                                     
    178       CASE ( -1 )                                ! esopa: test all possibility with control print 
    179                         CALL dyn_spg_exp( kt ) 
    180                         CALL prt_ctl( tab3d_1=ua, clinfo1=' spg0 - Ua: ', mask1=umask, & 
    181          &                            tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    182                         CALL dyn_spg_ts ( kt ) 
    183                         CALL prt_ctl( tab3d_1=ua, clinfo1=' spg1 - Ua: ', mask1=umask, & 
    184          &                           tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    185                         CALL dyn_spg_flt( kt, kindic ) 
    186                         CALL prt_ctl( tab3d_1=ua, clinfo1=' spg2 - Ua: ', mask1=umask, & 
    187          &                            tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    188168      END SELECT 
    189169      !                     
    190170      IF( l_trddyn )   THEN                      ! save the surface pressure gradient trends for further diagnostics 
    191          SELECT CASE ( nspg ) 
    192          CASE ( 0, 1 ) 
    193             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    194             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    195          CASE( 2 ) 
    196             z2dt = 2. * rdt 
    197             IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt 
    198             ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / z2dt - ztrdu(:,:,:) 
    199             ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:) 
    200          END SELECT 
     171         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     172         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    201173         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt ) 
    202174         ! 
     
    216188      !!                  ***  ROUTINE dyn_spg_init  *** 
    217189      !!                 
    218       !! ** Purpose :   Control the consistency between cpp options for  
     190      !! ** Purpose :   Control the consistency between namelist options for  
    219191      !!              surface pressure gradient schemes 
    220192      !!---------------------------------------------------------------------- 
    221       INTEGER ::   ioptio 
     193      INTEGER ::   ioptio, ios 
     194      ! 
     195      NAMELIST/namdyn_spg/ ln_dynspg_exp, ln_dynspg_ts, & 
     196      &                    ln_bt_fw, ln_bt_av, ln_bt_auto, & 
     197      &                    nn_baro, rn_bt_cmax, nn_bt_flt 
    222198      !!---------------------------------------------------------------------- 
    223199      ! 
    224200      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_init') 
    225201      ! 
    226       IF(lwp) THEN             ! Control print 
     202      REWIND( numnam_ref )              ! Namelist namdyn_spg in reference namelist : Free surface 
     203      READ  ( numnam_ref, namdyn_spg, IOSTAT = ios, ERR = 901) 
     204901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in reference namelist', lwp ) 
     205 
     206      REWIND( numnam_cfg )              ! Namelist namdyn_spg in configuration namelist : Free surface 
     207      READ  ( numnam_cfg, namdyn_spg, IOSTAT = ios, ERR = 902 ) 
     208902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist', lwp ) 
     209      IF(lwm) WRITE ( numond, namdyn_spg ) 
     210      ! 
     211      IF(lwp) THEN             ! Namelist print 
    227212         WRITE(numout,*) 
    228213         WRITE(numout,*) 'dyn_spg_init : choice of the surface pressure gradient scheme' 
    229214         WRITE(numout,*) '~~~~~~~~~~~' 
    230          WRITE(numout,*) '     Explicit free surface                  lk_dynspg_exp = ', lk_dynspg_exp 
    231          WRITE(numout,*) '     Free surface with time splitting       lk_dynspg_ts  = ', lk_dynspg_ts 
    232          WRITE(numout,*) '     Filtered free surface cst volume       lk_dynspg_flt = ', lk_dynspg_flt 
    233       ENDIF 
    234  
    235       IF( lk_dynspg_ts ) CALL dyn_spg_ts_init( nit000 ) 
    236       ! (do it now, to set nn_baro, used to allocate some arrays later on) 
    237       !                        ! allocate dyn_spg arrays 
    238       IF( lk_dynspg_ts ) THEN 
    239          IF( dynspg_oce_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_oce arrays') 
     215         WRITE(numout,*) '     Explicit free surface                  ln_dynspg_exp = ', ln_dynspg_exp 
     216         WRITE(numout,*) '     Free surface with time splitting       ln_dynspg_ts  = ', ln_dynspg_ts 
     217      ENDIF 
     218 
     219      IF( ln_dynspg_ts ) THEN 
     220         CALL dyn_spg_ts_init( nit000  ) ! do it first, to set nn_baro, used to allocate some arrays later on 
     221         !                               ! allocate dyn_spg arrays 
    240222         IF( dyn_spg_ts_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays') 
    241223         IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' ) 
     
    244226      !                        ! Control of surface pressure gradient scheme options 
    245227      ioptio = 0 
    246       IF(lk_dynspg_exp)   ioptio = ioptio + 1 
    247       IF(lk_dynspg_ts )   ioptio = ioptio + 1 
    248       IF(lk_dynspg_flt)   ioptio = ioptio + 1 
    249       ! 
    250       IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   & 
    251            &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 
    252       IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav )   & 
    253            &   CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' ) 
    254       ! 
    255       IF( lk_esopa     )   nspg = -1 
    256       IF( lk_dynspg_exp)   nspg =  0 
    257       IF( lk_dynspg_ts )   nspg =  1 
    258       IF( lk_dynspg_flt)   nspg =  2 
    259       ! 
    260       IF( lk_esopa     )   nspg = -1 
     228      IF(ln_dynspg_exp)   ioptio = ioptio + 1 
     229      IF(ln_dynspg_ts )   ioptio = ioptio + 1 
     230      ! 
     231      IF(  ioptio > 1  .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   & 
     232           &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme' ) 
     233      IF( ln_dynspg_ts .AND. ln_isfcav )   & 
     234           &   CALL ctl_stop( ' dynspg_ts not tested with ice shelf cavity ' ) 
     235      ! 
     236      IF( ln_dynspg_exp)   nspg =  0 
     237      IF( ln_dynspg_ts )   nspg =  1 
    261238      ! 
    262239      IF(lwp) THEN 
    263240         WRITE(numout,*) 
    264          IF( nspg == -1 )   WRITE(numout,*) '     ESOPA test All scheme used' 
    265241         IF( nspg ==  0 )   WRITE(numout,*) '     explicit free surface' 
    266242         IF( nspg ==  1 )   WRITE(numout,*) '     free surface with time splitting scheme' 
    267          IF( nspg ==  2 )   WRITE(numout,*) '     filtered free surface' 
    268       ENDIF 
    269  
    270 #if defined key_dynspg_flt || defined key_esopa 
    271       CALL solver_init( nit000 )   ! Elliptic solver initialisation 
    272 #endif 
    273  
    274       !                        ! Control of timestep choice 
    275       IF( lk_dynspg_ts .OR. lk_dynspg_exp ) THEN 
    276          IF( nn_cla == 1 )   CALL ctl_stop( 'Crossland advection not implemented for this free surface formulation' ) 
    277       ENDIF 
    278  
    279       !               ! Control of hydrostatic pressure choice 
    280       IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) THEN 
    281          CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 
    282243      ENDIF 
    283244      ! 
Note: See TracChangeset for help on using the changeset viewer.