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 6360 – NEMO

Changeset 6360


Ignore:
Timestamp:
2016-03-07T14:56:33+01:00 (8 years ago)
Author:
hliu
Message:

W/D updating

Location:
branches/2015/dev_r4826_NOC_WAD/NEMOGCM
Files:
20 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/CONFIG/AMM12/cpp_AMM12.fcm

    r4245 r6360  
    1  bld::tool::fppkeys  key_bdy key_tide key_dynspg_ts key_ldfslp  key_zdfgls  key_vvl key_diainstant key_mpp_mpi key_iomput 
     1 bld::tool::fppkeys  key_bdy key_dynspg_ts key_ldfslp  key_zdfgls  key_vvl key_mpp_mpi  
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/CONFIG/SHARED/namelist_ref

    r4773 r6360  
    11961196   rn_htrmax         =  200.0   ! max. depth of transition range 
    11971197/ 
     1198!----------------------------------------------------------------------- 
     1199&namwad  !       Wetting and Drying parameter setting  
     1200!----------------------------------------------------------------------- 
     1201  ln_wd         = .false.          !  with/out Wetting/Drying (T: on, F: off) 
     1202  rn_wdmin1     = 0.10             !  minimum water depth on dried land cells  
     1203  rn_wdmin2     = 0.01             !  tolerance of rn_wdmin1 (in iterations) 
     1204  rn_wdld       = 10.0             !  height measured from MSL where the W/D 
     1205                                   !  zone meets true land zone 
     1206  nn_wdit       = 5                !  maximum number of iteration for W/D 
     1207                                   !  limiter 
     1208/  
     1209 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/CONFIG/cfg.txt

    r4690 r6360  
    77C1D_PAPA OPA_SRC 
    88ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
    9 AMM12 OPA_SRC 
    109GYRE_BFM OPA_SRC TOP_SRC 
    1110ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    1211ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 
     12AMM12 OPA_SRC 
     13THACKER OPA_SRC 
     14IRS18WE OPA_SRC 
     15IRS18WD OPA_SRC 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90

    r4292 r6360  
    344344CONTAINS 
    345345   SUBROUTINE dia_hth( kt )         ! Empty routine 
     346      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    346347      WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt 
    347348   END SUBROUTINE dia_hth 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r5014 r6360  
    267267 
    268268   !!---------------------------------------------------------------------- 
    269    !! critical depths,limiters,and masks for  Wetting and Drying 
     269   !! critical depths,filters, limiters,and masks for  Wetting and Drying 
    270270   !! --------------------------------------------------------------------- 
    271271 
     272   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wduflt, wdvflt     !: u- and v- filter 
    272273   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wdmask             !: u- and v- limiter  
    273274 
     
    413414 
    414415      IF(ln_wd) & 
    415       ALLOCATE( wdmask(jpi,jpj), STAT=ierr(12) ) 
     416      ALLOCATE( wduflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj), STAT=ierr(12) ) 
    416417      ! 
    417418      dom_oce_alloc = MAXVAL(ierr) 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5441 r6360  
    718718                     &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    719719                     &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     720 
     721                  !!debugging for W/D 
     722                  !IF(pe3_out(ji,jj,jk) < 0._wp .AND. jk == 1) Then 
     723                  !!IF(ji+nimpp-1 == 21 .AND. jj+njmpp-1 == 283) THEN 
     724                  !  write(500+nproc,*), '---dom_vvl_interpol: negative value---' 
     725                  !  write(500+nproc,*), 'pout, ji,jj,jk,pe3_out,umask(ji,jj,jk)' 
     726                  !  write(500+nproc,*), pout, ji+nimpp-1, jj+njmpp-1, jk, pe3_out(ji,jj,jk), umask(ji,jj,jk) 
     727                  !  write(500+nproc,*), 'pe3_in, pe3_inp1, e3t_0, e3t_0p1' 
     728                  !  write(500+nproc,*),  pe3_in(ji,jj,jk), pe3_in(ji+1,jj,jk) 
     729                  !  write(500+nproc,*), e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) 
     730                  !  write(500+nproc,*), 'e12t,e12tp1,r1_e12u' 
     731                  !  write(500+nproc,*), e12t(ji,jj), e12t(ji+1,jj),r1_e12u(ji,jj) 
     732                  !  write(500+nproc,*), e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) 
     733                  !  write(500+nproc,*), e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) )  
     734                  !!  stop 
     735                  !END IF 
    720736               END DO 
    721737            END DO 
     
    809825      END SELECT 
    810826      ! 
     827 
     828      !!debugging for W/D 
     829      !   DO jk = 1, jpk 
     830      !      DO jj = 1, jpjm1 
     831      !         DO ji = 1, fs_jpim1   ! vector opt. 
     832      !            IF(pe3_out(ji,jj,jk) < 0._wp) Then 
     833      !              write(400+nproc,*), '---dom_vvl_interpol: negative value---' 
     834      !              write(400+nproc,*), 'pout, ji,jj,jk,pe3_out,umask(ji,jj,jk)' 
     835      !              write(400+nproc,*), pout, ji+nimpp-1, jj+njmpp-1, jk, pe3_out(ji,jj,jk), umask(ji,jj,jk) 
     836      !              write(400+nproc,*), 'pe3_in, pe3_inp1, e3t_0, e3t_0p1, e3u_0' 
     837      !              write(400+nproc,*),  pe3_in(ji,jj,jk), pe3_in(ji+1,jj,jk) 
     838      !              write(400+nproc,*), e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk), e3u_0(ji,jj,jk) 
     839      !              write(400+nproc,*), 'e12t,e12tp1,r1_e12u' 
     840      !              write(400+nproc,*), e12t(ji,jj), e12t(ji+1,jj),r1_e12u(ji,jj) 
     841      !              write(400+nproc,*), e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) 
     842      !              write(400+nproc,*), e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) )  
     843      !              !stop 
     844      !            END IF 
     845      !         END DO 
     846      !      END DO 
     847      !   END DO 
     848 
    811849 
    812850      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol') 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5466 r6360  
    12891289        DO jj = 1, jpj 
    12901290           DO ji = 1, jpi 
    1291               zenv(ji,jj) = MAX( zenv(ji,jj), -rn_wdld )    ! filt out land bathy data  
     1291              !zenv(ji,jj) = MAX( zenv(ji,jj), -rn_wdld )    ! filt out land bathy data  
    12921292           END DO 
    12931293        END DO 
     
    14901490      fse3uw(:,:,:) = e3uw_0  (:,:,:) 
    14911491      fse3vw(:,:,:) = e3vw_0  (:,:,:) 
     1492 
     1493        ! ! for tracing down the W/D stable issue 
     1494        ! ji=25-nimpp+1;jj=263-njmpp+1 
     1495        ! if (ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj )then 
     1496        ! write(800,'(2i6, 3f10.2)') ji, jj, bathy(ji,jj), bathy(ji+1, jj), bathy(ji, jj+1) 
     1497        ! write(800,*) 'jk,fsdeptij,depwij, de3wij,fsdeptip1j,depwip1j, de3wip1j,fsdeptijp1,depwijp1, de3wijp1' 
     1498        ! DO jk = 1, jpkm1 
     1499        !   write(800,'(i3,10e11.3)')  jk,fsdept(ji,jj,jk),fsdepw(ji,jj,jk),fsde3w(ji,jj,jk),& 
     1500        !                              &  fsdept(ji+1,jj,jk),fsdepw(ji+1,jj,jk),fsde3w(ji+1,jj,jk),& 
     1501        !                              &  fsdept(ji,jj+1,jk),fsdepw(ji,jj+1,jk),fsde3w(ji,jj+1,jk) 
     1502        ! 
     1503        ! END DO 
     1504        ! end if 
    14921505!! 
    14931506      ! HYBRID :  
     
    15301543      ENDIF 
    15311544      !  END DO 
    1532       IF(lwp) THEN                                  ! selected vertical profiles 
    1533          WRITE(numout,*) 
    1534          WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 
    1535          WRITE(numout,*) ' ~~~~~~  --------------------' 
    1536          WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
    1537          WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     & 
    1538             &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk ) 
    1539          DO jj = mj0(20), mj1(20) 
    1540             DO ji = mi0(20), mi1(20) 
    1541                WRITE(numout,*) 
    1542                WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj) 
    1543                WRITE(numout,*) ' ~~~~~~  --------------------' 
    1544                WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
    1545                WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     & 
    1546                   &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 
    1547             END DO 
    1548          END DO 
    1549          DO jj = mj0(74), mj1(74) 
    1550             DO ji = mi0(100), mi1(100) 
    1551                WRITE(numout,*) 
    1552                WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj) 
    1553                WRITE(numout,*) ' ~~~~~~  --------------------' 
    1554                WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
    1555                WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     & 
    1556                   &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 
    1557             END DO 
    1558          END DO 
    1559       ENDIF 
    1560  
     1545!      IF(lwp) THEN                                  ! selected vertical profiles 
     1546!         WRITE(numout,*) 
     1547!         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 
     1548!         WRITE(numout,*) ' ~~~~~~  --------------------' 
     1549!         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
     1550!         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     & 
     1551!            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk ) 
     1552!         DO jj = mj0(20), mj1(20) 
     1553!            DO ji = mi0(20), mi1(20) 
     1554!               WRITE(numout,*) 
     1555!               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj) 
     1556!               WRITE(numout,*) ' ~~~~~~  --------------------' 
     1557!               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
     1558!               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     & 
     1559!                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 
     1560!            END DO 
     1561!        END DO 
     1562!        DO jj = mj0(74), mj1(74) 
     1563!           DO ji = mi0(100), mi1(100) 
     1564!              WRITE(numout,*) 
     1565!              WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj) 
     1566!              WRITE(numout,*) ' ~~~~~~  --------------------' 
     1567!              WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')") 
     1568!              WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     & 
     1569!                 &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 
     1570!           END DO 
     1571!        END DO 
     1572!      ENDIF 
     1573! 
    15611574!================================================================================ 
    15621575! check the coordinate makes sense 
    15631576!================================================================================ 
    1564       IF(ln_wd) THEN 
    1565       DO ji = 1, jpi 
    1566          DO jj = 1, jpj 
    1567  
    1568             IF( bathy(ji,jj) > 0._wp) THEN 
    1569                DO jk = 1, mbathy(ji,jj) 
    1570                  ! check coordinate is monotonically increasing 
    1571                  IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
    1572                     WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    1573                     WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    1574                     WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 
    1575                     WRITE(numout,*) 'e3t',fse3t(ji,jj,:) 
    1576                     CALL ctl_stop( ctmp1 ) 
    1577                  ENDIF 
    1578                  ! and check it has never gone negative 
    1579                  IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 
    1580                     WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
    1581                     WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk 
    1582                     WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
    1583                     WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
    1584                     CALL ctl_stop( ctmp1 ) 
    1585                  ENDIF 
    1586                  ! and check it never exceeds the total depth 
    1587                  IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN 
    1588                     WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
    1589                     WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
    1590                     WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
    1591                     CALL ctl_stop( ctmp1 ) 
    1592                  ENDIF 
    1593                END DO 
    1594  
    1595                DO jk = 1, mbathy(ji,jj)-1 
    1596                  ! and check it never exceeds the total depth 
    1597                 IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN 
    1598                     WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
    1599                     WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
    1600                     WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
    1601                     CALL ctl_stop( ctmp1 ) 
    1602                  ENDIF 
    1603                END DO 
    1604  
    1605             ENDIF 
    1606  
    1607          END DO 
    1608       END DO 
    1609       END IF 
     1577 !     DO ji = 1, jpi 
     1578 !        DO jj = 1, jpj 
     1579 ! 
     1580 !           IF( bathy(ji,jj) > 0._wp) THEN 
     1581 !              DO jk = 1, mbathy(ji,jj) 
     1582 !                ! check coordinate is monotonically increasing 
     1583 !                IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
     1584 !                   WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1585 !                   WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1586 !                   WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 
     1587 !                   WRITE(numout,*) 'e3t',fse3t(ji,jj,:) 
     1588 !                   CALL ctl_stop( ctmp1 ) 
     1589 !                ENDIF 
     1590 !                ! and check it has never gone negative 
     1591 !                IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 
     1592 !                   WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
     1593 !                   WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1594 !                   WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
     1595 !                   WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
     1596 !                   CALL ctl_stop( ctmp1 ) 
     1597 !                ENDIF 
     1598 !                ! and check it never exceeds the total depth 
     1599 !                IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     1600 !                   WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1601 !                   WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1602 !                   WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
     1603 !                   CALL ctl_stop( ctmp1 ) 
     1604 !                ENDIF 
     1605 !              END DO 
     1606 ! 
     1607 !              DO jk = 1, mbathy(ji,jj)-1 
     1608 !                ! and check it never exceeds the total depth 
     1609 !               IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     1610 !                   WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1611 !                   WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1612 !                   WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
     1613 !                   CALL ctl_stop( ctmp1 ) 
     1614 !                ENDIF 
     1615 !              END DO 
     1616 ! 
     1617 !            ENDIF 
     1618 ! 
     1619 !         END DO 
     1620 !      END DO 
    16101621      ! 
    16111622      CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r4370 r6360  
    210210      ! 
    211211      DO jk = 1, jpk 
    212          tsn(:,:,jk,jp_tem) = (  ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) )   & 
    213             &                + 10. * ( 5000. - fsdept(:,:,jk) ) /5000.)  ) * tmask(:,:,jk) 
     212      !   tsn(:,:,jk,jp_tem) = (  ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) )   & 
     213      !      &                + 10. * ( 5000. - fsdept(:,:,jk) ) /5000.)  ) * tmask(:,:,jk) 
     214      !  for constant temperature testing 
     215         tsn(:,:,jk,jp_tem) = 10.0 
    214216         tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 
    215217      END DO 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90

    r3294 r6360  
    6969 
    7070 
     71        print*,'maximum value of bfru(v)a', maxval(bfrua), maxval(bfrva) 
    7172# if defined key_vectopt_loop 
    7273        DO jj = 1, 1 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r5221 r6360  
    5757   INTEGER , PUBLIC ::   nhpg  =  0   ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 
    5858 
     59   ! temporary debugging integer 
     60   INTEGER :: jidbg = 163, jjdbg = 179 
     61 
    5962   !! * Substitutions 
    6063#  include "domzgr_substitute.h90" 
     
    6568   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6669   !!---------------------------------------------------------------------- 
     70    
    6771CONTAINS 
    6872 
     
    369373      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    370374      !! 
    371       INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    372       REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp   ! temporary scalars 
    373       LOGICAL  ::   ll_tmp1, ll_tmp2           ! local logical variables 
     375      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices 
     376      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars 
     377      LOGICAL  ::   ll_tmp1, ll_tmp2, ll_tmp3            ! local logical variables 
    374378      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
    375       REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
     379      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy !W/D pressure filter 
    376380      !!---------------------------------------------------------------------- 
    377381      ! 
     
    384388         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used' 
    385389      ENDIF 
     390 
     391 
     392      ! debugging W/D  
     393      ! rhd(:,:,:) = 0.0_wp 
     394      ! end debugging w/D, must be removed later 
    386395 
    387396      ! Local constant initialization 
     
    395404        DO jj = 2, jpjm1 
    396405           DO ji = 2, jpim1  
    397              ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
    398              ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     406             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))  
     407             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2 
     408             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
    399409                                                       & rn_wdmin1 + rn_wdmin2 
    400410 
    401              IF(ll_tmp1) THEN 
     411             IF(ll_tmp1.AND.ll_tmp2) THEN 
    402412               zcpx(ji,jj) = 1.0_wp 
    403              ELSE IF(ll_tmp2) THEN 
     413               wduflt(ji,jj) = 1.0_wp 
     414             ELSE IF(ll_tmp3) THEN 
    404415               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
    405                zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     416               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / & 
    406417                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     418               zcpx(ji,jj) = MIN(zcpx(ji,jj), 1._wp) 
     419               wduflt(ji,jj) = 1.0_wp 
    407420             ELSE 
    408421               zcpx(ji,jj) = 0._wp 
     422               wduflt(ji,jj) = 0.0_wp 
    409423             END IF 
    410424       
    411              ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
    412              ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     425             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))  
     426             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2 
     427             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
    413428                                                       & rn_wdmin1 + rn_wdmin2 
    414429 
    415              IF(ll_tmp1) THEN 
     430             IF(ll_tmp1.AND.ll_tmp2) THEN 
    416431               zcpy(ji,jj) = 1.0_wp 
    417              ELSE IF(ll_tmp2) THEN 
     432               wdvflt(ji,jj) = 1.0_wp 
     433             ELSE IF(ll_tmp3) THEN 
    418434               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
    419                zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     435               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / & 
    420436                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     437               zcpy(ji,jj) = MIN(zcpy(ji,jj), 1._wp) 
     438               wdvflt(ji,jj) = 1.0_wp 
    421439             ELSE 
    422440               zcpy(ji,jj) = 0._wp 
     441               wdvflt(ji,jj) = 0.0_wp 
    423442             END IF 
    424443           END DO 
     
    426445        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    427446      ENDIF 
     447 
     448 
     449      jii=jidbg-nimpp+1;jjj=jjdbg-njmpp+1 
    428450 
    429451      ! Surface value 
     
    441463               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
    442464 
     465 
    443466            IF(ln_wd) THEN 
     467 
     468         if (jii>1 .and. jii <= jpim1 .and.  jjj>1 .and. jjj <= jpj )then 
     469         if(ji == jii .and. jj == jjj) then 
     470           write(700,*) 'bathy_ij, ip1j, ijp1', bathy(ji,jj), bathy(ji+1,jj), bathy(ji,jj+1) 
     471           write(600,*) 'fsde3wij,ip1j,ijp1, jk,sshnij,ip1j,ijp1', fsde3w(ji,jj,1), fsde3w(ji+1,jj,1),& 
     472           & fsde3w(ji,jj+1,1),1,sshn(ji,jj), sshn(ji+1,jj),sshn(ji,jj+1) 
     473           write(500,*) 'zhpi,zuap,zcpx,k', zhpi(ji,jj,1), zuap, zcpx(ji,jj), 1,& 
     474           &            ua(ji,jj,1),(zhpi(ji,jj,1) + zuap)*zcpx(ji,jj)  
     475           write(500,*) 'zhpj,zvap,zcpy,k', zhpj(ji,jj,1), zvap, zcpy(ji,jj), 1,& 
     476           &            va(ji,jj,1),(zhpj(ji,jj,1) + zvap)*zcpy(ji,jj) 
     477         end if 
     478         end if 
    444479              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
    445480              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     
    449484 
    450485            ! add to the general momentum trend 
     486            ! W/D testing line  
     487            IF(min(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1) THEN 
     488            !for debugging W/D, need to uncomment later 
    451489            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     490            END IF 
     491            IF(min(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1) THEN 
     492            !for debugging W/D, need to uncomment later 
    452493            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
     494            END IF 
     495 
     496            !end W/D debugging 
    453497         END DO 
    454498      END DO 
     
    472516 
    473517               IF(ln_wd) THEN 
     518         if (jii>1 .and. jii <= jpim1 .and.  jjj>1 .and. jjj <= jpj )then 
     519         if(ji == jii .and. jj == jjj) then 
     520           write(600,*) 'fsde3wij,ip1j,ijp1, jk,sshnij,ip1j,ijp1', fsde3w(ji,jj,jk), fsde3w(ji+1,jj,jk),& 
     521           &fsde3w(ji,jj+1,jk),jk,sshn(ji,jj), sshn(ji+1,jj),sshn(ji,jj+1) 
     522           write(500,*) 'zhpi,zuap,zcpx,k', zhpi(ji,jj,jk), zuap, zcpx(ji,jj), jk,& 
     523          &                                 ua(ji,jj,jk), (zhpi(ji,jj,jk) + zuap)*zcpx(ji,jj) 
     524           write(500,*) 'zhpj,zvap,zcpy,k', zhpj(ji,jj,jk), zvap, zcpy(ji,jj), jk,& 
     525          &                                 va(ji,jj,jk), (zhpj(ji,jj,jk) + zvap)*zcpy(ji,jj) 
     526         end if 
     527         end if 
    474528                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
    475529                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     
    479533 
    480534               ! add to the general momentum trend 
    481                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
    482                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
     535            ! W/D testing line  
     536            IF(min(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1) THEN 
     537            !for debugging W/D, need to uncomment later 
     538            ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     539            END IF 
     540            IF(min(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1) THEN 
     541            !for debugging W/D, need to uncomment later 
     542            va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
     543            END IF 
     544  
     545              !for debugging W/D, need to uncomment later 
     546              !ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     547               !va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
     548               !end W/D debugging 
    483549            END DO 
    484550         END DO 
    485551      END DO 
     552 
    486553      ! 
    487554      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     
    523590        DO jj = 2, jpjm1 
    524591           DO ji = 2, jpim1  
    525              ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     592             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 
     593                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 
     594                     &  > rn_wdmin1 + rn_wdmin2 
    526595             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
    527596                                                       & rn_wdmin1 + rn_wdmin2 
     
    537606             END IF 
    538607       
    539              ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     608             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
     609                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
     610                     &  > rn_wdmin1 + rn_wdmin2 
    540611             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
    541612                                                       & rn_wdmin1 + rn_wdmin2 
     
    819890        DO jj = 2, jpjm1 
    820891           DO ji = 2, jpim1  
    821              ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     892             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 
     893                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 
     894                     &  > rn_wdmin1 + rn_wdmin2 
    822895             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
    823896                                                       & rn_wdmin1 + rn_wdmin2 
     
    829902               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
    830903                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     904               zcpx(ji,jj) = MIN(zcpx(ji,jj), 1._wp) 
    831905             ELSE 
    832906               zcpx(ji,jj) = 0._wp 
    833907             END IF 
    834908       
    835              ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     909             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
     910                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
     911                     &  > rn_wdmin1 + rn_wdmin2 
    836912             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
    837913                                                       & rn_wdmin1 + rn_wdmin2 
     
    843919               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
    844920                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     921               zcpy(ji,jj) = MIN(zcpy(ji,jj), 1._wp) 
    845922             ELSE 
    846923               zcpy(ji,jj) = 0._wp 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90

    r3294 r6360  
    103103            END DO  
    104104         END DO 
     105 
     106         !! for tracing down the W/D stable issue 
     107         !ji=25-nimpp+1;jj=263-njmpp+1 
     108         !if (ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj .and. jk == 21 )then 
     109         !  write(300,*) 'zhke, advkeg-U', - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj), kt  
     110         !  write(300,*) 'zhke, advkeg-V', - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj), kt  
     111         !end if 
     112 
     113 
    105114!!gm idea to be tested  ==>>   is it faster on scalar computers ? 
    106115!         DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r5495 r6360  
    162162      REAL(wp), POINTER, DIMENSION(:,:) :: zhf 
    163163      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef. 
     164      REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1           ! Wetting/Dying velocity filter coef. 
    164165      !!---------------------------------------------------------------------- 
     166 
     167      ! tempoary debugging integer 
     168      INTEGER :: jidbg, jjdbg 
     169      jidbg = 106; jjdbg = 131 
    165170      ! 
    166171      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts') 
     
    173178      CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
    174179      CALL wrk_alloc( jpi, jpj, zhf ) 
    175       IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 
     180      IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
    176181      ! 
    177182      !                                         !* Local constant initialization 
     
    182187      zraur = 1._wp / rau0 
    183188      ! 
     189      ! 
    184190      IF( kt == nit000 .AND. neuler == 0 ) THEN         ! reciprocal of baroclinic time step 
    185191        z2dt_bf = rdt 
     
    228234            DO jj = 1, jpjm1 
    229235               DO ji = 1, jpim1 
    230                   zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                     & 
    231                         &          ht(ji  ,jj  ) + ht(ji+1,jj  )   )                   & 
    232                         &      / ( MAX( 1.0_wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    & 
     236                 !removed the inconsistancy in Jerome's version, hl 
     237                  zwz(ji,jj) =   ( ht(ji,jj+1)*tmask(ji,jj+1,1) + ht(ji+1,jj+1)*tmask(ji+1,jj+1,1) + & 
     238                        &          ht(ji,jj  )*tmask(ji,jj  ,1) + ht(ji+1,jj  )*tmask(ji+1,jj,  1) ) & 
     239                        &      / ( MAX( 1.0_wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +      & 
    233240                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) ) 
    234241                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj) 
     
    357364                &                                      + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
    358365                &                                      + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
     366                 IF(ji+nimpp-1 == jidbg .and. jj+njmpp-1 == jjdbg) THEN 
     367                   write(333,*) 'in spt_ts:  zutrdx, zutrdy calculations' 
     368                   write(333,*) 'U: ftne, zwy_ij,ftnwip1j, zwy_ip1j, ftseij,zwy_ijm1,ftswip1j, zwy_ip1jm1' 
     369                   write(333,*) ftne(ji,jj), zwy(ji,jj),ftnw(ji+1,jj), zwy(ji+1,jj)  
     370                   write(333,*) ftse(ji,jj), zwy(ji,jj-1),ftsw(ji+1,jj),zwy(ji+1,jj-1) 
     371                   write(333,*) 'V: ftsw, zwx_im1jp1,ftseijp1, zwx_ijp1, ftnwij,zwx_im1j,ftneij, zwx_ij' 
     372                   write(333,*) ftsw(ji,jj+1),zwx(ji-1,jj+1),ftse(ji,jj+1),  zwx(ji  ,jj+1) 
     373                   write(333,*)  ftnw(ji,jj), zwx(ji-1,jj), ftne(ji,jj), zwx(ji,jj) 
     374                 END IF 
     375 
    359376            END DO 
    360377         END DO 
     
    375392                IF(ll_tmp1) THEN 
    376393                  zcpx(ji,jj) = 1.0_wp 
     394                  wduflt1(ji,jj) = 1.0_wp 
    377395                ELSE IF(ll_tmp2) THEN 
    378396                   ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
    379397                  zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
    380398                              &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     399                  wduflt1(ji,jj) = 1.0_wp 
    381400                ELSE 
    382                    zcpx(ji,jj) = 0._wp 
     401                  zcpx(ji,jj) = 0._wp 
     402                  wduflt1(ji,jj) = 0.0_wp 
    383403                END IF 
    384404 
     405                !for w/d debugging, delete it when finished 
     406                !   zcpx(ji,jj) = 0._wp 
     407                ! end debugging part 
     408                    
    385409                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
    386410                        & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
     
    390414                IF(ll_tmp1) THEN 
    391415                   zcpy(ji,jj) = 1.0_wp 
     416                   wdvflt1(ji,jj) = 1.0_wp 
    392417                ELSE IF(ll_tmp2) THEN 
    393418                   ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
    394419                  zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
    395420                              &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     421                   wdvflt1(ji,jj) = 1.0_wp 
    396422                ELSE 
    397423                  zcpy(ji,jj) = 0._wp 
     424                  wdvflt1(ji,jj) = 0.0_wp 
    398425                END IF 
     426 
     427                !for w/d debugging, delete it when finished 
     428                !   zcpy(ji,jj) = 0._wp 
     429                ! end debugging part 
    399430             END DO 
    400431           END DO 
     
    405436           DO jj = 2, jpjm1 
    406437              DO ji = 2, jpim1 
    407                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) ) / & 
    408                                &                 e1u(ji,jj) * zcpx(ji,jj) 
    409                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) ) / & 
    410                                &                 e2v(ji,jj) * zcpy(ji,jj) 
     438                 IF(min(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1) THEN 
     439                 zu_trd(ji,jj) = (zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) ) / & 
     440                               &                 e1u(ji,jj)) * zcpx(ji,jj) 
     441                 END IF 
     442                 IF(min(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1) THEN 
     443                 zv_trd(ji,jj) = (zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) ) / & 
     444                               &                 e2v(ji,jj)) * zcpy(ji,jj) 
     445                  END IF 
     446 
     447                 IF(ji+nimpp-1 == jidbg .and. jj+njmpp-1 == jjdbg) THEN 
     448                   write(333,*) 'in spt_ts:  zutrdx, zutrdy, hpgx, hpgy, zcpx, zcpy' 
     449                   write(333,*) zu_trd(ji,jj), zv_trd(ji,jj), & 
     450                             &  -grav*(sshn(ji+1,jj)-sshn(ji,jj))/e1u(ji,jj)*zcpx(ji,jj), & 
     451                             &  -grav*(sshn(ji,jj+1)-sshn(ji,jj))/e2v(ji,jj)*zcpy(ji,jj), & 
     452                             &  zcpx(ji,jj), zcpy(ji,jj)  
     453                   write(334,*) 'in spt_ts:  sshn_ij, sshn_ip1j, sshn_ijp1, bathy_ij, bathy_ip1j, bathy_ijp1' 
     454                   write(334,*) sshn(ji,jj), sshn(ji+1,jj), sshn(ji,jj+1), bathy(ji,jj), bathy(ji+1,jj), bathy(ji,jj+1)  
     455                 END IF 
    411456              END DO 
    412457           END DO 
     
    415460              DO ji = fs_2, fs_jpim1   ! vector opt. 
    416461                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj) 
    417                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj) 
     462                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj)  
    418463              END DO 
    419464           END DO 
     
    422467      ENDIF 
    423468 
     469      IF(ln_wd) THEN 
     470      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend 
     471         DO ji = fs_2, fs_jpim1 
     472             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1) * wduflt1(ji,jj) 
     473             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) * wdvflt1(ji,jj) 
     474          END DO 
     475      END DO 
     476      ELSE 
    424477      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend 
    425478         DO ji = fs_2, fs_jpim1 
     
    428481          END DO 
    429482      END DO 
     483      END IF 
    430484      ! 
    431485      !                                         ! Add bottom stress contribution from baroclinic velocities: 
     
    489543            END DO 
    490544         ENDIF 
     545      ENDIF 
     546 
     547      IF ( ln_wd ) THEN                    ! Apply wetting/drying velocity filters  
     548            DO jj = 2, jpjm1 
     549               DO ji = fs_2, fs_jpim1   ! vector opt. 
     550                  zu_frc(ji,jj) = zu_frc(ji,jj) * wduflt1(ji,jj) 
     551                  zv_frc(ji,jj) = zv_frc(ji,jj) * wdvflt1(ji,jj) 
     552               END DO 
     553            END DO 
    491554      ENDIF 
    492555      !                                   !* Right-Hand-Side of the barotropic ssh equation 
     
    717780          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
    718781 
     782         IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
     783           DO jj = 2, jpjm1 
     784              DO ji = 2, jpim1 
     785                 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))& 
     786                        & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj)) & 
     787                        &  > rn_wdmin1 + rn_wdmin2 
     788                 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
     789                                                           & rn_wdmin1 + rn_wdmin2 
     790                 IF(ll_tmp1) THEN 
     791                   zcpx(ji,jj) = 1.0_wp 
     792                   wduflt1(ji,jj) = 1.0_wp 
     793                 ELSE IF(ll_tmp2) THEN 
     794                    ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
     795                   zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 
     796                               & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 
     797                   wduflt1(ji,jj) = 1.0_wp 
     798                 ELSE 
     799                    zcpx(ji,jj) = 0._wp 
     800                    wduflt1(ji,jj) = 0.0_wp 
     801                 END IF 
     802 
     803                !for w/d debugging, delete it when finished 
     804                !   zcpx(ji,jj) = 0._wp 
     805                ! end debugging part 
     806 
     807                 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))& 
     808                        & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1)) & 
     809                        &  > rn_wdmin1 + rn_wdmin2 
     810                 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
     811                                                           & rn_wdmin1 + rn_wdmin2 
     812                 IF(ll_tmp1) THEN 
     813                    zcpy(ji,jj) = 1.0_wp 
     814                    wdvflt1(ji,jj) = 1.0_wp 
     815                 ELSE IF(ll_tmp2) THEN 
     816                    ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
     817                   zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 
     818                               & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 
     819                    wdvflt1(ji,jj) = 1.0_wp 
     820                 ELSE 
     821                   zcpy(ji,jj) = 0._wp 
     822                    wdvflt1(ji,jj) = 0.0_wp 
     823                 END IF 
     824 
     825                !for w/d debugging, delete it when finished 
     826                !   zcpy(ji,jj) = 0._wp 
     827                ! end debugging part 
     828              END DO 
     829            END DO 
     830            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     831         ENDIF 
    719832         ! 
    720833         ! Compute associated depths at U and V points: 
     
    801914         ENDIF 
    802915         ! 
     916 
     917 
    803918         ! Add bottom stresses: 
    804919         IF(ln_wd) THEN 
     
    809924           zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 
    810925         END IF 
     926 
     927         ! filter out all dynamic forcing on land cells 
     928         IF(ln_wd) THEN 
     929           zu_trd(:,:) = zu_trd(:,:) * wduflt1(:,:) 
     930           zv_trd(:,:) = zv_trd(:,:) * wduflt1(:,:) 
     931         END IF 
    811932         ! 
    812933         ! Surface pressure trend: 
    813          IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
    814            DO jj = 2, jpjm1 
    815               DO ji = 2, jpim1 
    816                  ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))& 
    817                         & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj)) & 
    818                         &  > rn_wdmin1 + rn_wdmin2 
    819                  ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
    820                                                            & rn_wdmin1 + rn_wdmin2 
    821                  IF(ll_tmp1) THEN 
    822                    zcpx(ji,jj) = 1.0_wp 
    823                  ELSE IF(ll_tmp2) THEN 
    824                     ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
    825                    zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 
    826                                & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 
    827                  ELSE 
    828                     zcpx(ji,jj) = 0._wp 
    829                  END IF 
    830  
    831                  ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))& 
    832                         & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1)) & 
    833                         &  > rn_wdmin1 + rn_wdmin2 
    834                  ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
    835                                                            & rn_wdmin1 + rn_wdmin2 
    836                  IF(ll_tmp1) THEN 
    837                     zcpy(ji,jj) = 1.0_wp 
    838                  ELSE IF(ll_tmp2) THEN 
    839                     ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
    840                    zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 
    841                                & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 
    842                  ELSE 
    843                    zcpy(ji,jj) = 0._wp 
    844                  END IF 
    845               END DO 
    846             END DO 
    847             CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    848          ENDIF 
    849934 
    850935         IF(ln_wd) THEN 
     
    854939                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
    855940                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
    856                  zwx(ji,jj) = zu_spg * zcpx(ji,jj) 
     941                 ! for W/D debugging 
     942                 !zwx(ji,jj) = zu_spg * zcpx(ji,jj) * 0._wp 
     943                 !zwy(ji,jj) = zv_spg * zcpy(ji,jj) * 0._wp 
     944                 IF(ji+nimpp-1 == jidbg .and. jj+njmpp-1 == jjdbg) THEN 
     945                   write(444,*) 'in spt_ts_intg:  zu_spg, zv_spg, zcpx, zcpy' 
     946                   write(444,*) zu_spg, zv_spg, zcpx(ji,jj), zcpy(ji,jj) 
     947                   write(445,*) 'in spt_ts_intg:  sshn_ij, sshn_ip1j, sshn_ijp1, bathy_ij, bathy_ip1j, bathy_ijp1' 
     948                   write(445,*) zsshp2_e(ji,jj), zsshp2_e(ji+1,jj), zsshp2_e(ji,jj+1), & 
     949                    &           bathy(ji,jj), bathy(ji+1,jj), bathy(ji,jj+1)  
     950                 END IF 
     951                 zwx(ji,jj) = zu_spg * zcpx(ji,jj)  
    857952                 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 
    858953              END DO 
     
    879974                            &                                 + zu_trd(ji,jj)   & 
    880975                            &                                 + zu_frc(ji,jj) ) & 
    881                             &   ) * umask(ji,jj,1) 
     976                            &   * wduflt(ji,jj)) * umask(ji,jj,1) 
    882977 
    883978                  va_e(ji,jj) = (                                zvn_e(ji,jj)   & 
     
    885980                            &                                 + zv_trd(ji,jj)   & 
    886981                            &                                 + zv_frc(ji,jj) ) & 
    887                             &   ) * vmask(ji,jj,1) 
     982                            &   * wdvflt(ji,jj)) * vmask(ji,jj,1) 
    888983               END DO 
    889984            END DO 
     
    9081003                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   & 
    9091004                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) & 
    910                             &   ) * zhura 
     1005                            &   * wduflt(ji,jj)) * zhura 
    9111006 
    9121007                  va_e(ji,jj) = (                hv_e(ji,jj)  *  zvn_e(ji,jj)   & 
     
    9141009                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   & 
    9151010                            &               +      hv(ji,jj)  * zv_frc(ji,jj) ) & 
    916                             &   ) * zhvra 
     1011                            &   * wdvflt(ji,jj) ) * zhvra 
    9171012               END DO 
    9181013            END DO 
     
    10601155      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
    10611156      CALL wrk_dealloc( jpi, jpj, zhf ) 
    1062       IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 
     1157      IF(ln_wd) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
    10631158      ! 
    10641159      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
     
    11031198              CASE( 1 )  ! Boxcar, width = nn_baro 
    11041199                 DO jn = 1, 3*nn_baro 
    1105                     za1 = ABS(float(jn-jic))/float(nn_baro) 
     1200                    za1 = ABS(real(jn-jic,wp))/real(nn_baro,wp) 
    11061201                    IF (za1 < 0.5_wp) THEN 
    11071202                      zwgt1(jn) = 1._wp 
     
    11121207              CASE( 2 )  ! Boxcar, width = 2 * nn_baro 
    11131208                 DO jn = 1, 3*nn_baro 
    1114                     za1 = ABS(float(jn-jic))/float(nn_baro) 
     1209                    za1 = ABS(real(jn-jic,wp))/real(nn_baro,wp) 
    11151210                    IF (za1 < 1._wp) THEN 
    11161211                      zwgt1(jn) = 1._wp 
     
    12501345      IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 
    12511346 
    1252       rdtbt = rdt / FLOAT(nn_baro) 
     1347      rdtbt = rdt / real(nn_baro,wp) 
    12531348      zcmax = zcmax * rdtbt 
    12541349                                                        ! Print results 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r4624 r6360  
    645645            CALL lbc_lnk( zwz, 'F', 1. ) 
    646646        CASE ( 4 )                                                ! total (relative + planetary vorticity) 
     647         !ji=2-nimpp+1;jj=229-njmpp+1 
     648         !if ( ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj .and. jk == 1)then 
     649         !  write(400+nproc,*) 'in dynvor, rotn(ji,jj,1), ff(ji,jj), fse3f(ji,jj,1), ze3f(ji,jj,1), kt' 
     650         !  write(400+nproc,*) rotn(ji,jj,jk), ff(ji,jj), fse3f(ji,jj,1), ze3f(ji,jj,jk), kt 
     651         !end if 
    647652            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 
    648653         CASE ( 5 )                                                ! total (coriolis + metric) 
     
    680685            END DO 
    681686         END DO 
     687 
     688         !ji=2-nimpp+1;jj=229-njmpp+1 
     689         !if ( ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj .and. jk == 1)then 
     690         !  write(400+nproc,*) 'in dynvor, een, zwzij, zwzim1j, zwzijm1, zwzim1jm1, rotnij, kt' 
     691         !  write(400+nproc,*) zwz(ji,jj),zwz(ji-1,jj),zwz(ji,jj-1),zwz(ji-1,jj-1), rotn(ji,jj,1), kt 
     692         !end if 
     693 
    682694         DO jj = 2, jpjm1 
    683695            DO ji = fs_2, fs_jpim1   ! vector opt. 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

    r3294 r6360  
    102102      END DO 
    103103 
     104         !ji=25-nimpp+1;jj=263-njmpp+1;jk=21 
     105         !if (ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj )then 
     106         !  !write(400,*) 'advw-U', ua(ji,jj,jk) - ztrdu(ji,jj,jk) 
     107         !  !write(400,*) 'advw-V', va(ji,jj,jk) - ztrdv(ji,jj,jk)  
     108         !  write(400,*) 'advw-U1,zwuw,ukm1, uk', ua(ji,jj,jk), kt, zwuw(ji,jj,jk), un(ji,jj,jk-1), un(ji,jj,jk) 
     109         !  write(400,*) 'advw-V1,zwvw,vkm1, vk', va(ji,jj,jk), kt, zwvw(ji,jj,jk), vn(ji,jj,jk-1), vn(ji,jj,jk) 
     110         !end if 
     111 
    104112      DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points 
    105113         DO jj = 2, jpjm1 
     
    114122         END DO   
    115123      END DO 
     124 
     125         !ji=25-nimpp+1;jj=263-njmpp+1;jk=21 
     126         !if (ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj )then 
     127         !  !write(400,*) 'advw-U', ua(ji,jj,jk) - ztrdu(ji,jj,jk) 
     128         !  !write(400,*) 'advw-V', va(ji,jj,jk) - ztrdv(ji,jj,jk)  
     129         !  write(400,*) 'advw-U2', ua(ji,jj,jk), kt 
     130         !  write(400,*) 'advw-V2', va(ji,jj,jk), kt 
     131         !end if 
    116132 
    117133      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r4328 r6360  
    132132         DO jn = 1, jpts 
    133133            DO jk = 1, jpkm1 
    134                tsn(:,:,jk,jn) = tsa(:,:,jk,jn)     
     134               !tsn(:,:,jk,jn) = tsa(:,:,jk,jn)     
    135135            END DO 
    136136         END DO 
     
    318318 
    319319                   ze3t_f = 1.e0 / ze3t_f 
    320                    ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered 
    321                    ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta 
     320                   !ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered 
     321                   !ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta 
    322322                   ! 
    323323                   IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only) 
    324324                      ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d ) 
    325                       pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average 
     325                      !pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average 
    326326                   ENDIF 
    327327               END DO 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90

    r4605 r6360  
    597597   END SUBROUTINE trd_icp_init 
    598598   SUBROUTINE trd_dwr( kt )          ! Empty routine 
     599      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    599600      WRITE(*,*) 'trd_dwr: You should not have seen this print! error ?', kt 
    600601   END SUBROUTINE trd_dwr 
    601602   SUBROUTINE trd_twr( kt )          ! Empty routine 
     603      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    602604      WRITE(*,*) 'trd_twr: You should not have seen this print! error ?', kt 
    603605   END SUBROUTINE trd_twr 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r4624 r6360  
    153153         END DO 
    154154 
     155 
    155156         ! 
    156157         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition 
    157158         ! 
     159         !print*, 'MAXVALUE of bfru(v)a in zdfbfr----', maxval(bfrua), maxval(bfrva) 
    158160         IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        & 
    159161            &                       tab2d_2=bfrva, clinfo2=       ' v: ', mask2=vmask,ovlap=1 ) 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

    r4624 r6360  
    248248CONTAINS 
    249249   SUBROUTINE zdf_ddm( kt )           ! Dummy routine 
     250      INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step 
    250251      WRITE(*,*) 'zdf_ddm: You should not have seen this print! error?', kt 
    251252   END SUBROUTINE zdf_ddm 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r5135 r6360  
    262262      REWIND( numnam_ref )              ! Namelist namwad in reference namelist : Parameters for Wetting/Drying 
    263263      READ  ( numnam_ref, namwad, IOSTAT = ios, ERR = 905) 
    264 905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', lwp ) 
     264905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', .TRUE. ) 
    265265 
    266266      REWIND( numnam_cfg )              ! Namelist namwad in configuration namelist : Parameters for Wetting/Drying 
    267267      READ  ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906) 
    268 906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', lwp ) 
     268906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', .TRUE. ) 
    269269 
    270270 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/step.F90

    r4760 r6360  
    3131   !!---------------------------------------------------------------------- 
    3232   USE step_oce         ! time stepping definition modules 
     33   USE, intrinsic :: ieee_arithmetic 
    3334 
    3435   IMPLICIT NONE 
     
    7071      !!              -8- Outputs and diagnostics 
    7172      !!---------------------------------------------------------------------- 
    72       INTEGER ::   jk       ! dummy loop indice 
     73      INTEGER ::   ji,jj, jk       ! dummy loop indice 
    7374      INTEGER ::   indic    ! error indicator if < 0 
    7475      INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt) 
     76 
     77      !! temporary debug info 
     78      INTEGER ::   jidbg, jjdbg, jkdbg 
     79 
     80      jidbg = 163; jjdbg = 179 
    7581      !! --------------------------------------------------------------------- 
    7682 
     
    180186          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! subtract Neptune velocities (simplified) 
    181187          IF( lk_bdy           )  CALL bdy_dyn3d_dmp( kstp )   ! bdy damping trends 
    182                                   CALL dyn_adv      ( kstp )   ! advection (vector or flux form) 
     188 
     189          ! tracing down trend trouble 
     190                                  ji=jidbg-nimpp+1;jj=jjdbg-njmpp+1 
     191                                  if (ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj )then 
     192                                  write(300+nproc,*) 'adv1',kstp,ua(ji,jj,1),va(ji,jj,1),umask(ji,jj,1) 
     193                                  end if 
     194 
     195!                                  IF(count(ieee_is_nan(ua(:,:,:)))+count(ieee_is_nan(va(:,:,:)) > 0 )) THEN 
     196!                                  write(300+nproc, *) 'before calling dyn_adv in time step:', kstp 
     197!                                  DO jk = 1, jpkm1; DO jj = 1, jpjm1; DO ji = 1, jpim1 
     198!                                    IF (ieee_is_nan(ua(ji,jj,jk)+va(ji,jj,jk))) THEN 
     199!   write(300+nproc,'(5i5, 4f10.2)') ji,jj,ji+nimpp-1, jj+njmpp-1,jk,ua(ji,jj,jk), va(ji,jj,jk), bathy(ji,jj), sshn(ji,jj) 
     200!                                      STOP 
     201!                                    ENDIF 
     202!                                  END DO; END DO; END DO 
     203!                                  END IF 
     204 
     205                                  CALL dyn_adv      ( kstp )   !  advection (vector or flux form) 
     206 
     207                                  ji=jidbg-nimpp+1;jj=jjdbg-njmpp+1 
     208                                  if ( ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj )then 
     209                                     write(300+nproc,*) 'adv2',kstp,ua(ji,jj,1),va(ji,jj,1) 
     210                                  endif  
     211 
     212!                                  IF(count(ieee_is_nan(ua(:,:,:)))+count(ieee_is_nan(va(:,:,:)) > 0 )) THEN 
     213!                                  write(300+nproc, *) 'after calling dyn_adv in time step:', kstp 
     214!                                  DO jk = 1, jpkm1; DO jj = 1, jpjm1; DO ji = 1, jpim1 
     215!                                    IF (ieee_is_nan(ua(ji,jj,jk)+va(ji,jj,jk))) THEN 
     216!   write(300+nproc,'(5i5, 4f10.2)') ji,jj,ji+nimpp-1, jj+njmpp-1,jk,ua(ji,jj,jk), va(ji,jj,jk), bathy(ji,jj), sshn(ji,jj) 
     217!                                      STOP 
     218!                                    ENDIF 
     219!                                  END DO; END DO; END DO 
     220!                                  END IF 
     221 
    183222                                  CALL dyn_vor      ( kstp )   ! vorticity term including Coriolis 
     223                                   
     224                                  ji=jidbg-nimpp+1;jj=jjdbg-njmpp+1 
     225                                  if ( ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj )then 
     226                                    write(300+nproc,*) 'cor',ua(ji,jj,1),va(ji,jj,1) 
     227                                  endif  
     228 
     229!                                  IF(count(ieee_is_nan(ua(:,:,:)))+count(ieee_is_nan(va(:,:,:)) > 0 )) THEN 
     230!                                  write(300+nproc, *) 'after calling dyn_vor in time step:', kstp 
     231!                                  DO jk = 1, jpkm1; DO jj = 1, jpjm1; DO ji = 1, jpim1 
     232!                                    IF (ieee_is_nan(ua(ji,jj,jk)+va(ji,jj,jk))) THEN 
     233!   write(300+nproc,'(5i5, 4f10.2)') ji,jj,ji+nimpp-1, jj+njmpp-1,jk,ua(ji,jj,jk), va(ji,jj,jk), bathy(ji,jj), sshn(ji,jj) 
     234!                                      STOP 
     235!                                    ENDIF 
     236!                                  END DO; END DO; END DO 
     237!                                  END IF 
     238 
    184239                                  CALL dyn_ldf      ( kstp )   ! lateral mixing 
     240!                                  IF(count(ieee_is_nan(ua(:,:,:)))+count(ieee_is_nan(va(:,:,:)) > 0 )) THEN 
     241!                                  write(300+nproc, *) 'after calling dyn_ldf in time step:', kstp 
     242!                                  DO jk = 1, jpkm1; DO jj = 1, jpjm1; DO ji = 1, jpim1 
     243!                                    IF (ieee_is_nan(ua(ji,jj,jk)+va(ji,jj,jk))) THEN 
     244!   write(300+nproc,'(5i5, 4f10.2)') ji,jj,ji+nimpp-1, jj+njmpp-1,jk,ua(ji,jj,jk), va(ji,jj,jk), bathy(ji,jj), sshn(ji,jj) 
     245!                                      STOP 
     246!                                    ENDIF 
     247!                                  END DO; END DO; END DO 
     248!                                  END IF 
     249 
    185250          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! add Neptune velocities (simplified) 
    186251#if defined key_agrif 
    187252          IF(.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn        ! momentum sponge 
    188253#endif 
     254                                  ji=jidbg-nimpp+1;jj=jjdbg-njmpp+1 
     255                                  if ( ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj )then 
     256                                     write(300+nproc,*) 'ldf',ua(ji,jj,1),va(ji,jj,1) 
     257                                  endif  
     258 
    189259                                  CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure 
     260!                                  IF(count(ieee_is_nan(ua(:,:,:)))+count(ieee_is_nan(va(:,:,:)) > 0 )) THEN 
     261!                                  write(300+nproc, *) 'after calling dyn_hpg in time step:', kstp 
     262!                                  DO jk = 1, jpkm1; DO jj = 1, jpjm1; DO ji = 1, jpim1 
     263!                                    IF (ieee_is_nan(ua(ji,jj,jk)+va(ji,jj,jk))) THEN 
     264!   write(300+nproc,'(5i5, 4f10.2)') ji,jj,ji+nimpp-1, jj+njmpp-1,jk,ua(ji,jj,jk), va(ji,jj,jk), bathy(ji,jj), sshn(ji,jj) 
     265!                                      STOP 
     266!                                    ENDIF 
     267!                                  END DO; END DO; END DO 
     268!                                  END IF 
     269 
     270                                  ji=jidbg-nimpp+1;jj=jjdbg-njmpp+1 
     271                                  if (ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj )then 
     272                                    write(300+nproc,*) 'hpg',ua(ji,jj,1),va(ji,jj,1) 
     273                                  endif  
     274 
    190275                                  CALL dyn_spg( kstp, indic )  ! surface pressure gradient 
     276                                  if ( ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj )then 
     277                                   write(300+nproc,*) 'spg',ua(ji,jj,1),va(ji,jj,1) 
     278                                  endif 
     279 
     280!                                  IF(count(ieee_is_nan(ua(:,:,:)))+count(ieee_is_nan(va(:,:,:)) > 0 )) THEN 
     281!                                  write(300+nproc, *) 'after calling dyn_spg in time step:', kstp 
     282!                                  DO jk = 1, jpkm1; DO jj = 1, jpjm1; DO ji = 1, jpim1 
     283!                                    IF (ieee_is_nan(ua(ji,jj,jk)+va(ji,jj,jk))) THEN 
     284!   write(300+nproc,'(5i5, 4f10.2)') ji,jj,ji+nimpp-1, jj+njmpp-1,jk,ua(ji,jj,jk), va(ji,jj,jk), bathy(ji,jj), sshn(ji,jj) 
     285!                                      STOP 
     286!                                    ENDIF 
     287!                                  END DO; END DO; END DO 
     288!                                  END IF 
    191289 
    192290                                  ua_sv(:,:,:) = ua(:,:,:)     ! Save trends (barotropic trend has been fully updated at this stage) 
     
    258356      ENDIF 
    259357 
     358                                  IF(count(ieee_is_nan(ua(:,:,:)))+count(ieee_is_nan(va(:,:,:)) > 0 )) THEN 
     359                                  write(300+nproc, *) 'after calling tracer part in time step:', kstp 
     360                                  DO jk = 1, jpkm1; DO jj = 1, jpjm1; DO ji = 1, jpim1 
     361                                    IF (ieee_is_nan(ua(ji,jj,jk)+va(ji,jj,jk))) THEN 
     362   write(300+nproc,'(5i5, 4f10.2)') ji,jj,ji+nimpp-1, jj+njmpp-1,jk,ua(ji,jj,jk), va(ji,jj,jk), bathy(ji,jj), sshn(ji,jj) 
     363                                      STOP 
     364                                    ENDIF 
     365                                  END DO; END DO; END DO 
     366                                  END IF 
    260367      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    261368      ! Dynamics                                    (tsa used as workspace) 
     
    271378                               rotn(:,:,:)  = rotb(:,:,:)  
    272379 
     380                                !ji=25-nimpp+1;jj=263-njmpp+1 
     381                                !if (ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj )then 
     382                                !  write(300+nproc,*) 'dyn_bfr1',ua(ji,jj,1),ua(ji,jj,21),ji,jpi 
     383                                !endif  
    273384                               CALL dyn_bfr( kstp )         ! bottom friction 
     385                               ! if (ji>1 .and. ji <= jpim1 .and.  jj>1 .and. jj <= jpj )then 
     386                               !   write(300+nproc,*) 'dyn_bfr2',ua(ji,jj,1),ua(ji,jj,21),ji,jpi 
     387                               ! endif  
    274388                               CALL dyn_zdf( kstp )         ! vertical diffusion 
    275389      ELSE 
     
    294408                               CALL dyn_spg( kstp, indic )  ! surface pressure gradient 
    295409      ENDIF 
     410                               IF(ln_wd) THEN 
     411                                 DO jk = 1, jpkm1 
     412                                   ua(:,:,jk) = ua(:,:,jk) * wduflt(:,:) 
     413                                   va(:,:,jk) = va(:,:,jk) * wdvflt(:,:) 
     414                                 END DO 
     415                               END IF 
     416 
    296417                               CALL dyn_nxt( kstp )         ! lateral velocity at next time step 
    297418 
Note: See TracChangeset for help on using the changeset viewer.