Ignore:
Timestamp:
2014-12-02T16:13:53+01:00 (6 years ago)
Author:
rblod
Message:

dev_r4642_WavesWG : stokes-coriolis from sbcwave + time splitting

Location:
branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r4624 r4960  
    3939   USE wrk_nemo        ! Memory Allocation 
    4040   USE timing          ! Timing     
     41   USE sbcwave 
    4142   USE sbcapr          ! surface boundary condition: atmospheric pressure 
    4243   USE dynadv, ONLY: ln_dynadv_vec 
     
    412413         zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:) 
    413414      ENDIF   
     415 
     416      ! Add Stokes Coriolis if defined 
     417      IF ( ln_stcor ) THEN 
     418         DO jj = 1, jpjm1 
     419            DO ji = 1, fs_jpim1   ! vector opt. 
     420 
     421               zy1 = ff(ji  ,jj-1) * ( vsd2d(ji  ,jj-1) + vsd2d(ji+1,jj-1) ) 
     422               zy2 = ff(ji  ,jj  ) * ( vsd2d(ji  ,jj  ) + vsd2d(ji+1,jj  ) ) 
     423               zx1 = ff(ji-1,jj  ) * ( usd2d(ji-1,jj  ) + usd2d(ji-1,jj+1) ) 
     424               zx2 = ff(ji  ,jj  ) * ( usd2d(ji  ,jj  ) + usd2d(ji  ,jj+1) ) 
     425 
     426               zu_frc(ji,jj) = zu_frc(ji,jj) + 0.25 * ( zy1 + zy2 ) * hur(ji,jj) 
     427               zv_frc(ji,jj) = zv_frc(ji,jj) - 0.25 * ( zx1 + zx2 ) * hvr(ji,jj) 
     428           END DO 
     429         END DO 
     430      ENDIF 
     431 
    414432      ! 
    415433      IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing 
  • branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90

    r4644 r4960  
    1818   USE wrk_nemo         ! Memory Allocation 
    1919   USE sbcmod           ! Access to ln_stcor (sbc_oce) and wave parameters (sbc_wave) 
    20    USE sbcwave_ecmwf    ! Wave module 
     20   USE sbcwave          ! Wave module 
     21   USE timing          ! Timing 
    2122 
    2223   IMPLICIT NONE 
    23  
    24    REAL(wp) :: rn_deptmaxstcor = 150.0_wp ! maximum depth [m] to be affected by Stokes-Coriolis effect 
    25    !PRIVATE 
     24   PRIVATE 
    2625 
    2726   !! * Routine accessibility 
    28    PUBLIC dyn_stcor           ! routine called by step.F90 
    29  
    30    !! * Shared module variables 
    31    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ustc, vstc  ! Stokes-Coriolis u and v 
    32  
    33    !! * Module variables 
     27   PUBLIC dyn_stcor      ! routine called by step.F90 
    3428 
    3529   !! * Substitutions 
     
    3731#  include "domzgr_substitute.h90" 
    3832   !!---------------------------------------------------------------------- 
    39    !!   OPA 9.0 , implemented by Bedford Institute of Oceanography 
     33   !! NEMO/OPA 3.6 , NEMO Consortium (2010) 
     34   !! $Id:$ 
     35   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4036   !!---------------------------------------------------------------------- 
    41  
    4237 CONTAINS 
    43  
    44    INTEGER FUNCTION dynstcor_alloc() 
    45       !!---------------------------------------------------------------------- 
    46       !!                    ***  ROUTINE dynstcor_alloc  *** 
    47       !! 
    48       !! History :  1.0  !  2012-10  (Oyvind Breivik) 
    49       !!---------------------------------------------------------------------- 
    50       ALLOCATE( ustc(jpi,jpj,jpk) , vstc(jpi,jpj,jpk) , STAT=dynstcor_alloc ) 
    51          ! 
    52       IF( dynstcor_alloc /= 0 )   CALL ctl_warn('dynstcor_alloc: array allocate failed.') 
    53    END FUNCTION dynstcor_alloc 
    5438 
    5539   SUBROUTINE dyn_stcor( kt ) 
     
    5943      !! ** Purpose:  Add Stokes-Coriolis forcing to horizontal momentum equation. 
    6044      !! 
    61       !! ** History:  0.1  !   2012-10  oyvind.breivik@ecmwf.int 
     45      !! ** History:  3.4  !   2012-10  (O. Breivik)   Initial version 
     46      !!              3.6  !   2014-10  (R. Benshila)  
    6247      !!---------------------------------------------------------------------- 
    6348      INTEGER, INTENT( in ) ::   kt  ! ocean time-step index 
    6449      !! 
    6550      INTEGER  ::  ji, jj, jk        ! dummy loop indices 
    66       REAL(wp) ::  ztransp, zsp0, zk, zfac 
    67       REAL(wp) ::  zus, zvs 
     51      REAL(wp) :: zx1, zx2, zy1, zy2 
     52      !!---------------------------------------------------------------------- 
    6853      ! 
    69  
    70       ! Allocation at first time step. 
    71  
    72       IF ( kt == nit000 ) THEN 
    73          IF( dynstcor_alloc() /= 0 ) CALL ctl_stop('dyn_stcor: array allocate failed.') 
     54      IF( nn_timing == 1 )  CALL timing_start('dyn_stcor') 
     55      ! 
     56      IF( kt == nit000 ) THEN 
     57         IF(lwp) WRITE(numout,*) 
     58         IF(lwp) WRITE(numout,*) 'dyn_stcor : time stepping' 
     59         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    7460      ENDIF 
    7561 
    76       !!---------------------------------------------------------------------- 
    7762      ! 
    78       ! Update velocity tendencies ua, va by adding the Stokes-Coriolis velocities ustc, vstc 
     63      ! Update velocity tendencies ua, va by adding the Stokes-Coriolis velocities  
    7964      ! 
    80       DO jk = 1, jpk 
    81          DO jj = 1, jpj 
    82             DO ji = 1, jpi 
    83                ! Skip deep levels where the Stokes-Coriolis effect is negligible 
    84                IF (fsdept(ji,jj,jk) <= rn_deptmaxstcor) THEN 
    85                   ! Stokes transport speed estimated from Hs and Tmean 
    86                   ztransp = 2.0_wp*rpi*swh_wavepar(ji,jj)**2.0_wp/(16.0_wp*MAX(mwp_wavepar(ji,jj),0.0000001_wp)) 
     65      DO jk = 1, jpkm1 
     66         DO jj = 2, jpjm1 
     67            DO ji = 2, jpim1 
     68               zy1 = ff(ji  ,jj-1) * ( vsd3d(ji  ,jj-1,jk) + vsd3d(ji+1,jj-1,jk) ) 
     69               zy2 = ff(ji  ,jj  ) * ( vsd3d(ji  ,jj  ,jk) + vsd3d(ji+1,jj  ,jk) ) 
     70               zx1 = ff(ji-1,jj  ) * ( usd3d(ji-1,jj  ,jk) + usd3d(ji-1,jj+1,jk) ) 
     71               zx2 = ff(ji  ,jj  ) * ( usd3d(ji  ,jj  ,jk) + usd3d(ji  ,jj+1,jk) ) 
    8772 
    88                   ! Stokes surface speed 
    89                   zsp0 = SQRT(ust_wavepar(ji,jj)**2 + vst_wavepar(ji,jj)**2) 
    90  
    91                   ! Wavenumber scale 
    92                   zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp) 
    93  
    94                   ! Depth attenuation 
    95                   zfac = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk)) 
    96  
    97                   ! The Stokes-Coriolis forcing 
    98                   zus =  ff(ji,jj)*vst_wavepar(ji,jj)*zfac 
    99                   zvs = -ff(ji,jj)*ust_wavepar(ji,jj)*zfac 
    100  
    101                   ! Store arrays of tendencies for diagnostics output 
    102                   ! This may be removed later for efficiency 
    103                   ustc(ji,jj,jk) = zus 
    104                   vstc(ji,jj,jk) = zvs 
    105  
    106                   ua(ji,jj,jk) = ua(ji,jj,jk) + zus * umask(ji,jj,jk) 
    107                   va(ji,jj,jk) = va(ji,jj,jk) + zvs * vmask(ji,jj,jk) 
    108                ENDIF 
     73               ua(ji,jj,jk) = ua(ji,jj,jk) + 0.25 * ( zy1 + zy2 ) *umask(ji,jj,jk) 
     74               va(ji,jj,jk) = va(ji,jj,jk) - 0.25 * ( zx1 + zx2 ) *vmask(ji,jj,jk) 
    10975            ENDDO 
    11076         ENDDO 
    11177      ENDDO 
    11278      ! 
     79      IF( nn_timing == 1 )  CALL timing_stop('dynst_cor') 
     80      ! 
    11381   END SUBROUTINE dyn_stcor 
    11482 
    115  
    11683END MODULE dynstcor 
  • branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r4644 r4960  
    324324      ! Read wave parameters if Stokes-Coriolis forcing OR wave parameters for TKE are needed 
    325325      ! Also read it if drag coefficient to ensure that we use the neutral 10m from WAM. 
    326       IF (ln_stcor .OR. ln_wavetke .OR. ln_tauoc .OR. ln_cdgw) THEN 
     326      IF (ln_wavetke .OR. ln_tauoc .OR. ln_cdgw) THEN 
    327327         CALL sbc_wavepar( kt ) 
    328328      ENDIF 
  • branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r4644 r4960  
    3131   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
    3232   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
    33    REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum  
     33   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum  
    3434   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d  
    3535 
     
    144144            END DO 
    145145         END DO 
     146         ! 
     147         CALL lbc_lnk( usd2d(:,:), 'U', -1. ) 
     148         CALL lbc_lnk( vsd2d(:,:), 'V', -1. ) 
    146149 
    147           !Computation of the 3d Stokes Drift 
    148           DO jk = 1, jpk 
    149              DO jj = 1, jpj-1 
    150                 DO ji = 1, jpi-1 
    151                    usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk)))) 
    152                    vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk)))) 
    153                 END DO 
    154              END DO 
    155              usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) ) 
    156              vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,jpj,jk)) ) 
    157           END DO 
     150         !Computation of the 3d Stokes Drift 
     151         DO jk = 1, jpk 
     152            DO jj = 1, jpj-1 
     153               DO ji = 1, jpi-1 
     154                  usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk)))) 
     155                  vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk)))) 
     156               END DO 
     157            END DO 
     158            usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) ) 
     159            vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,jpj,jk)) ) 
     160         END DO 
     161         ! 
     162         CALL lbc_lnk( usd3d(:,:,:), 'U', -1. ) 
     163         CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. ) 
    158164 
    159165          CALL wrk_alloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 
  • branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step.F90

    r4644 r4960  
    184184                                  CALL dyn_ldf      ( kstp )   ! lateral mixing 
    185185          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! add Neptune velocities (simplified) 
     186          IF( ln_stcor    )       CALL dyn_stcor    ( kstp )   ! Stokes-Coriolis forcing (ln_stcor set in SBC/sbc_oce.F90) 
    186187#if defined key_agrif 
    187188          IF(.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn        ! momentum sponge 
Note: See TracChangeset for help on using the changeset viewer.