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 7422 for branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 – NEMO

Ignore:
Timestamp:
2016-12-01T18:17:41+01:00 (7 years ago)
Author:
gm
Message:

#1805 dev_INGV_UKMO_2016: improve Stokes drift (including dynspg_ts , Stokes-Coriolis force , and GLS surface roughness

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r7351 r7422  
    1515   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility 
    1616   !!             3.7  ! 2015-11  (J. Chanut) free surface simplification 
     17   !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence 
    1718   !!--------------------------------------------------------------------- 
    1819 
     
    3839   USE sbctide         ! tides 
    3940   USE updtide         ! tide potential 
     41   USE sbcwave         ! surface wave 
    4042   ! 
    4143   USE in_out_manager  ! I/O manager 
     
    168170      CALL wrk_alloc( jpi,jpj,   zsshu_a, zsshv_a                  ) 
    169171      CALL wrk_alloc( jpi,jpj,   zhf ) 
    170       IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
    171       ! 
    172       zmdi=1.e+20                               !  missing data indicator for masking 
     172      IF( ln_wd )   CALL wrk_alloc( jpi,jpj,   zcpx, zcpy, wduflt1, wdvflt1 ) 
     173      ! 
    173174      !                                         !* Local constant initialization 
    174       z1_12 = 1._wp / 12._wp  
     175      zmdi = 1.e+20                                   !  missing data indicator for masking 
     176      ! 
     177      z1_12 = 1._wp / 12._wp                          ! constants 
    175178      z1_8  = 0.125_wp                                    
    176179      z1_4  = 0.25_wp 
    177180      z1_2  = 0.5_wp      
    178181      zraur = 1._wp / rau0 
    179       !                                            ! reciprocal of baroclinic time step  
     182      !                                               ! baroclinic time step  
    180183      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt 
    181184      ELSE                                        ;   z2dt_bf = 2.0_wp * rdt 
    182185      ENDIF 
    183       z1_2dt_b = 1.0_wp / z2dt_bf  
    184       ! 
    185       ll_init     = ln_bt_av                       ! if no time averaging, then no specific restart  
     186      z1_2dt_b = 1.0_wp / z2dt_bf                     ! reciprocal of baroclinic time step  
     187      ! 
     188      ll_init     = ln_bt_av                          ! if no time averaging, then no specific restart  
    186189      ll_fw_start = .FALSE. 
    187       !                                            ! time offset in steps for bdy data update 
     190      !                                               ! time offset in steps for bdy data update 
    188191      IF( .NOT.ln_bt_fw ) THEN   ;   noffset = - nn_baro 
    189192      ELSE                       ;   noffset =   0  
     
    255258            zwz(:,:) = 0._wp 
    256259            zhf(:,:) = 0._wp 
    257             IF ( .not. ln_sco ) THEN 
     260            IF( .not. ln_sco ) THEN 
    258261 
    259262!!gm  agree the JC comment  : this should be done in a much clear way 
     
    324327         END DO 
    325328      END DO 
     329       
     330!!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum.... 
     331!!gm  Is it correct to do so ?   I think so... 
     332       
     333       
    326334      !                                   !* barotropic Coriolis trends (vorticity scheme dependent) 
    327335      !                                   ! -------------------------------------------------------- 
     
    373381      !                                   ! ---------------------------------------------------- 
    374382      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient 
    375         IF( ln_wd ) THEN                        ! Calculating and applying W/D gravity filters 
    376           wduflt1(:,:) = 1.0_wp 
    377           wdvflt1(:,:) = 1.0_wp 
    378           DO jj = 2, jpjm1 
    379              DO ji = 2, jpim1 
    380                 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   & 
     383         IF( ln_wd ) THEN                        ! Calculating and applying W/D gravity filters 
     384            wduflt1(:,:) = 1.0_wp 
     385            wdvflt1(:,:) = 1.0_wp 
     386            DO jj = 2, jpjm1 
     387               DO ji = 2, jpim1 
     388                  ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   & 
    381389                        & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj))   & 
    382390                        &  > rn_wdmin1 + rn_wdmin2 
    383                 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   & 
     391                  ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   & 
    384392                        &                                   + rn_wdmin1 + rn_wdmin2 
    385                 IF(ll_tmp1) THEN 
    386                   zcpx(ji,jj)    = 1.0_wp 
    387                 ELSEIF(ll_tmp2) THEN 
    388                    ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen here 
    389                   zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 
     393                  IF(ll_tmp1) THEN 
     394                     zcpx(ji,jj)    = 1.0_wp 
     395                  ELSEIF(ll_tmp2) THEN 
     396                     ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen here 
     397                     zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 
    390398                        &          /(sshn(ji+1,jj) - sshn(ji,jj))) 
    391                 ELSE 
    392                   zcpx(ji,jj)    = 0._wp 
    393                   wduflt1(ji,jj) = 0.0_wp 
    394                 END IF 
    395  
    396                 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
     399                  ELSE 
     400                     zcpx(ji,jj)    = 0._wp 
     401                     wduflt1(ji,jj) = 0.0_wp 
     402                  END IF 
     403 
     404                  ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
    397405                        & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1))   & 
    398406                        &  > rn_wdmin1 + rn_wdmin2 
    399                 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
     407                  ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
    400408                        &                                   + rn_wdmin1 + rn_wdmin2 
    401                 IF(ll_tmp1) THEN 
    402                    zcpy(ji,jj)    = 1.0_wp 
    403                 ELSEIF(ll_tmp2) THEN 
    404                    ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen here 
    405                   zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 
     409                  IF(ll_tmp1) THEN 
     410                     zcpy(ji,jj)    = 1.0_wp 
     411                  ELSEIF(ll_tmp2) THEN 
     412                     ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen here 
     413                     zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 
    406414                        &          /(sshn(ji,jj+1) - sshn(ji,jj))) 
    407                 ELSE 
    408                   zcpy(ji,jj)    = 0._wp 
    409                   wdvflt1(ji,jj) = 0.0_wp 
    410                 ENDIF 
    411  
    412              END DO 
    413            END DO 
    414  
    415            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    416  
    417            DO jj = 2, jpjm1 
    418               DO ji = 2, jpim1 
    419                  zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
    420                         &                        * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj) 
    421                  zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
    422                         &                        * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj) 
    423               END DO 
    424            END DO 
    425  
    426          ELSE 
    427  
    428            DO jj = 2, jpjm1 
    429               DO ji = fs_2, fs_jpim1   ! vector opt. 
    430                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj) 
    431                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj)  
    432               END DO 
    433            END DO 
    434         ENDIF 
    435  
    436       ENDIF 
     415                  ELSE 
     416                     zcpy(ji,jj)    = 0._wp 
     417                     wdvflt1(ji,jj) = 0.0_wp 
     418                  ENDIF 
     419 
     420               END DO 
     421            END DO 
     422            ! 
     423            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     424            ! 
     425            DO jj = 2, jpjm1 
     426               DO ji = 2, jpim1 
     427                  zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
     428                         &                        * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj) 
     429                  zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
     430                         &                        * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj) 
     431               END DO 
     432            END DO 
     433            ! 
     434         ELSE        ! no Wet & Drying 
     435            ! 
     436            DO jj = 2, jpjm1 
     437               DO ji = fs_2, fs_jpim1   ! vector opt. 
     438                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj) 
     439                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj)  
     440               END DO 
     441            END DO 
     442         ENDIF 
     443         ! 
     444      ENDIF          ! end non linear free surface case 
    437445 
    438446      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend 
     
    473481      END IF 
    474482      ! 
     483!!gm   TOP stress only !!!   this should be with a test on ISF use or not 
    475484      !                                         ! Add top stress contribution from baroclinic velocities:       
    476       IF (ln_bt_fw) THEN 
     485      IF( ln_bt_fw ) THEN 
    477486         DO jj = 2, jpjm1 
    478487            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    538547                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     ) 
    539548      ENDIF 
     549      ! 
     550      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary 
     551         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 
     552      ENDIF 
     553      ! 
    540554#if defined key_asminc 
    541555      !                                         ! Include the IAU weighted SSH increment 
Note: See TracChangeset for help on using the changeset viewer.