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 4380 for branches/2013 – NEMO

Changeset 4380 for branches/2013


Ignore:
Timestamp:
2014-01-29T14:54:00+01:00 (10 years ago)
Author:
hliu
Message:

wetting and drying: some bugs removed from dynspg_ts.F90

Location:
branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r4375 r4380  
    372372            h_oce     = gdepw_0(jpk) 
    373373         ELSE                                         ! bump centered in the basin 
     374            !IF(lwp) WRITE(numout,*) 
     375            !IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump' 
     376            !ii_bump = jpidta / 2                           ! i-index of the bump center 
     377            !ij_bump = jpjdta / 2                           ! j-index of the bump center 
     378            !r_bump  = 50000._wp                            ! bump radius (meters)        
     379            !h_bump  =  2700._wp                            ! bump height (meters) 
     380            !h_oce   = gdepw_0(jpk)                         ! background ocean depth (meters) 
     381            !IF(lwp) WRITE(numout,*) '            bump characteristics: ' 
     382            !IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump 
     383            !IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters' 
     384            !IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index' 
     385            !IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters' 
     386            !!                                         
     387            !DO jj = 1, jpjdta                              ! zdta : 
     388            !   DO ji = 1, jpidta 
     389            !      zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump 
     390            !      zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump 
     391            !      zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) ) 
     392            !   END DO 
     393            !END DO 
     394 
    374395            IF(lwp) WRITE(numout,*) 
    375             IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump' 
    376             ii_bump = jpidta / 2                           ! i-index of the bump center 
    377             ij_bump = jpjdta / 2                           ! j-index of the bump center 
    378             r_bump  = 50000._wp                            ! bump radius (meters)        
    379             h_bump  =  2700._wp                            ! bump height (meters) 
     396            IF(lwp) WRITE(numout,*) '         bathymetry field: Thacker''s parabolic basin' 
     397            ii_bump = jpidta / 2 + 1                       ! i-index of the basin center 
     398            ij_bump = jpjdta / 2 + 1                       ! j-index of the basin center 
     399            r_bump  = 430620._wp                           ! basin radius (meters)        
     400            h_bump  =     50._wp                           ! basin depth (meters) 
    380401            h_oce   = gdepw_0(jpk)                         ! background ocean depth (meters) 
    381             IF(lwp) WRITE(numout,*) '            bump characteristics: ' 
    382             IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump 
    383             IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters' 
    384             IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index' 
     402            IF(lwp) WRITE(numout,*) '           basin characteristics: ' 
     403            IF(lwp) WRITE(numout,*) '              basin center (i,j)   = ', ii_bump, ii_bump 
     404            IF(lwp) WRITE(numout,*) '              basin depth          = ', h_bump , ' meters' 
     405            IF(lwp) WRITE(numout,*) '              basin radius         = ', r_bump , ' index' 
    385406            IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters' 
    386             !                                         
     407                                                     
    387408            DO jj = 1, jpjdta                              ! zdta : 
    388409               DO ji = 1, jpidta 
    389                   zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump 
    390                   zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump 
    391                   zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) ) 
     410                  zi = FLOAT( ji - ii_bump ) * ppe1_m  
     411                  zj = FLOAT( jj - ij_bump ) * ppe2_m  
     412                  zdta(ji,jj) = h_bump * ( 1._wp - ( zi*zi + zj*zj ) / (r_bump * r_bump) ) 
    392413               END DO 
    393414            END DO 
     415            !IF(lwp) WRITE(numout,*) 
     416            !IF(lwp) WRITE(numout,*) '         bathymetry field: Thacker''s parabolic channel' 
     417            !ii_bump = jpidta / 2                           ! i-index of the bump center 
     418            !r_bump  = 81000._wp                            ! bump radius (meters)        
     419            !h_bump  =    20._wp                            ! bump height (meters) 
     420            !h_oce   = gdepw_0(jpk)                         ! background ocean depth (meters) 
     421            !IF(lwp) WRITE(numout,*) '            channel characteristics: ' 
     422            !IF(lwp) WRITE(numout,*) '               channel center (i,j)   = ', ii_bump, ii_bump 
     423            !IF(lwp) WRITE(numout,*) '               channel depth          = ', h_bump , ' meters' 
     424            !IF(lwp) WRITE(numout,*) '               channel radius         = ', r_bump , ' index' 
     425            !IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters' 
     426            !                                         
     427            !DO jj = 1, jpjdta                              ! zdta : 
     428            !   DO ji = 1, jpidta 
     429            !      zi = FLOAT( ji - ii_bump ) * ppe1_m  
     430            !      zdta(ji,jj) = h_bump * ( 1._wp - (zi/r_bump) ** 2) - 10._wp 
     431            !   END DO 
     432            !END DO 
    394433            !                                              ! idta : 
    395434            IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk 
     
    14181457                  IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) THEN 
    14191458                    mbathy(ji,jj) = MAX( 2, jk ) 
    1420                   ELSE IF(scobot(ji,jj) <= rn_landele) THEN 
     1459                  ELSE IF(scobot(ji,jj) <= - rn_landele) THEN 
    14211460                    mbathy(ji,jj) = 0 
    14221461                  ELSE 
  • branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r4376 r4380  
    127127      REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum 
    128128 
    129       REAL(wp), POINTER, DIMENSION(:,:) :: zwadflt 
     129      REAL(wp), POINTER, DIMENSION(:,:) :: zwadfltu, zwadfltv 
    130130      !!---------------------------------------------------------------------- 
    131131      ! 
     
    136136      CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 
    137137 
    138       IF(ln_wad) CALL wrk_alloc( jpi, jpj, zwadflt) 
     138      IF(ln_wad) CALL wrk_alloc( jpi, jpj, zwadfltu, zwadfltv) 
    139139      ! 
    140140      IF( kt == nit000 ) THEN             !* initialisation 
     
    185185      ENDIF 
    186186 
    187       IF(ln_wad) zwadflt(:,:) = 1._wp 
     187      IF(ln_wad) THEN 
     188        zwadfltu(:,:) = 1._wp 
     189        zwadfltv(:,:) = 1._wp 
     190      END IF 
    188191 
    189192      ! ----------------------------------------------------------------------------- 
     
    412415         !                                                !* after ssh_e 
    413416         !                                                !  ----------- 
    414          DO jj = 2, jpjm1                                 ! Horizontal divergence of barotropic transports 
    415             DO ji = fs_2, fs_jpim1   ! vector opt. 
    416                zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     & 
    417                   &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     & 
    418                   &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     & 
    419                   &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) ) 
    420             END DO 
    421          END DO 
     417         IF(ln_wad) THEN 
     418            DO jj = 2, jpjm1                                 ! Horizontal divergence of barotropic transports 
     419               DO ji = fs_2, fs_jpim1   ! vector opt. 
     420                  zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj) * zwadfltu(ji,   jj  )   & 
     421                     &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj) * zwadfltu(ji-1, jj  )   & 
     422                     &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  ) * zwadfltv(ji,   jj  )   & 
     423                     &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1) * zwadfltv(ji,   jj-1) ) & 
     424                     &         / ( e1t(ji,jj) * e2t(ji,jj) ) 
     425               END DO 
     426            END DO 
     427         ELSE 
     428            DO jj = 2, jpjm1                                 ! Horizontal divergence of barotropic transports 
     429               DO ji = fs_2, fs_jpim1   ! vector opt. 
     430                  zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     & 
     431                     &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     & 
     432                     &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     & 
     433                     &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) & 
     434                     &         / ( e1t(ji,jj) * e2t(ji,jj) ) 
     435               END DO 
     436            END DO 
     437         END IF 
    422438         ! 
    423439#if defined key_obc 
     
    433449#endif 
    434450         ! 
     451         DO jj = 2, jpjm1                                      ! leap-frog on ssh_e 
     452            DO ji = fs_2, fs_jpim1   ! vector opt. 
     453               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * & 
     454                             &   ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)  
     455            END DO 
     456         END DO 
     457 
     458         !! generate W/D filter 
    435459         IF(ln_wad) THEN 
    436             DO jj = 2, jpjm1                                      ! leap-frog on ssh_e 
     460            DO jj = 2, jpjm1 
    437461               DO ji = fs_2, fs_jpim1   ! vector opt. 
    438                   ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * & 
    439                                 &   ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)  
    440                   IF(ssha_e(ji,jj) <= rn_wadmin) THEN 
    441                     zwadflt(ji,   jj  ) = 0._wp 
    442                     zwadflt(ji-1, jj  ) = 0._wp 
    443                     zwadflt(ji,   jj-1) = 0._wp 
    444                     zwadflt(ji-1, jj-1) = 0._wp 
     462                  IF(ssha_e(ji,jj) + bathy(ji,jj) <= rn_wadmin) THEN 
     463                    zwadfltu(ji,   jj  ) = 0._wp 
     464                    zwadfltu(ji-1, jj  ) = 0._wp 
     465                    zwadfltv(ji,   jj  ) = 0._wp 
     466                    zwadfltv(ji,   jj-1) = 0._wp 
    445467                  END IF 
    446468               END DO 
    447469            END DO 
    448          ELSE 
    449             DO jj = 2, jpjm1                                      ! leap-frog on ssh_e 
    450                DO ji = fs_2, fs_jpim1   ! vector opt. 
    451                   ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * & 
    452                                 &   ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)  
    453                END DO 
    454             END DO 
     470            CALL lbc_lnk(zwadfltu, 'U', 1._wp)  
     471            CALL lbc_lnk(zwadfltv, 'V', 1._wp)  
    455472         END IF 
     473 
    456474 
    457475         !                                                !* after barotropic velocities (vorticity scheme dependent) 
     
    622640#endif 
    623641 
     642         IF(ln_wad) THEN 
     643            DO jj = 2, jpjm1 
     644               DO ji = fs_2, fs_jpim1   ! vector opt. 
     645                 ua_e(ji,jj) = ua_e(ji,jj) * zwadfltu(ji,jj) 
     646                 va_e(ji,jj) = va_e(ji,jj) * zwadfltv(ji,jj) 
     647               END DO 
     648            END DO 
     649         END IF 
     650 
    624651         ! 
    625652         CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries  
     
    693720      IF(ln_wad) THEN 
    694721         DO jk=1,jpkm1 
    695             ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b * zwadflt(:,:) 
    696             va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b * zwadflt(:,:) 
     722            ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b * zwadfltu(:,:) 
     723            va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b * zwadfltv(:,:) 
    697724         END DO 
    698725      ELSE 
     
    714741      CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 
    715742 
    716       IF(ln_wad) CALL wrk_dealloc( jpi, jpj, zwadflt) 
     743      IF(ln_wad) CALL wrk_dealloc( jpi, jpj, zwadfltu, zwadfltv) 
    717744      ! 
    718745      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
Note: See TracChangeset for help on using the changeset viewer.