Changeset 5216


Ignore:
Timestamp:
2015-04-15T20:54:28+02:00 (6 years ago)
Author:
mathiot
Message:

ISF branch: minor changes in zpshde, sbcisf + bug correction in domzgr (e3wu) only if ice shelf

Location:
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5200 r5216  
    969969      !! 
    970970      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
    971       INTEGER  ::   ik, it          ! temporary integers 
     971      INTEGER  ::   ik, it, ikb, ikt ! temporary integers 
    972972      LOGICAL  ::   ll_print         ! Allow  control print for debugging 
    973973      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
     
    10931093         e3vw_0(:,:,jk) = e3w_1d(jk) 
    10941094      END DO 
     1095 
    10951096      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
    10961097         DO jj = 1, jpjm1 
     
    11051106      IF ( ln_isfcav ) THEN 
    11061107      ! (ISF) define e3uw (adapted for 2 cells in the water column) 
    1107       ! Need to test if the modification of only mikt and mbkt levels is enough 
    1108          DO jk = 2,jpk                           
    1109             DO jj = 1, jpjm1  
    1110                DO ji = 1, fs_jpim1   ! vector opt.  
    1111                   e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) & 
    1112                     &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) ) 
    1113                   e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) & 
    1114                     &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) ) 
    1115                END DO  
    1116             END DO  
     1108         DO jj = 2, jpjm1  
     1109            DO ji = 2, fs_jpim1   ! vector opt.  
     1110               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 
     1111               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 
     1112               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) & 
     1113                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) ) 
     1114               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 
     1115               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 
     1116               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) & 
     1117                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) ) 
     1118            END DO 
    11171119         END DO 
    11181120      END IF 
    1119        
     1121 
    11201122      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
    11211123      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    11221124      ! 
     1125 
    11231126      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    11241127         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
     
    12631266         risfdep(:,:) = 0. ; misfdep(:,:) = 1 
    12641267      END WHERE 
     1268 
     1269      ! remove very shallow ice shelf (less than ~ 10m if 75L) 
     1270      IF ( cp_cfg .NE. "isomip" ) THEN 
     1271         WHERE (misfdep(:,:) <= 10 .AND. misfdep(:,:) .GT. 1) 
     1272            misfdep = 0; risfdep = 0.0_wp; 
     1273            mbathy  = 0; bathy   = 0.0_wp; 
     1274         END WHERE 
     1275         WHERE (bathy(:,:) <= 30.0 .AND. gphit < -60) 
     1276            misfdep = 0; risfdep = 0.0_wp; 
     1277            mbathy  = 0; bathy   = 0.0_wp; 
     1278         END WHERE 
     1279      END IF 
    12651280  
    12661281! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 
     
    17871802               ENDIF  
    17881803            !       ... on ik / ik-1  
    1789                e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
     1804               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
    17901805               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    17911806! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r5204 r5216  
    585585 
    586586      !! compute ustar 
    587          zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:)) ) 
     587         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 
    588588 
    589589      !! Compute gammats 
     
    596596      !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 
    597597      !! compute ustar 
    598          zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:)) ) 
     598         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 
    599599 
    600600      !! compute Pr and Sc number (can be improved) 
     
    678678 
    679679      REAL(wp) :: ze3, zhk 
    680       REAL(wp), DIMENSION(:,:), POINTER :: zikt 
     680      REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl 
    681681 
    682682      INTEGER :: ji,jj,jk 
    683683      INTEGER :: ikt,ikb 
    684       INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb 
    685  
    686       CALL wrk_alloc( jpi,jpj, mkt, mkb  ) 
    687       CALL wrk_alloc( jpi,jpj, zikt ) 
     684 
     685      CALL wrk_alloc( jpi,jpj, zhisf_tbl) 
    688686 
    689687      ! get first and last level of tbl 
    690       mkt(:,:) = misfkt(:,:) 
    691       mkb(:,:) = misfkb(:,:) 
    692688 
    693689      varout(:,:)=0._wp 
    694       DO jj = 2,jpj 
    695          DO ji = 2,jpi 
    696             ikt = mkt(ji,jj) 
    697             ikb = mkb(ji,jj) 
     690      IF (cptin == 'U') THEN 
     691         DO jj = 1,jpj 
     692            DO ji = 1,jpi 
     693               ikt = miku(ji,jj) ; ikb = miku(ji,jj) 
     694            ! thickness of boundary layer at least the top level thickness 
     695               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt)) 
     696 
     697            ! determine the deepest level influenced by the boundary layer 
     698            ! test on tmask useless ????? 
     699               DO jk = ikt, mbku(ji,jj) 
     700                  IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 
     701               END DO 
     702               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
    698703 
    699704            ! level fully include in the ice shelf boundary layer 
    700             DO jk = ikt, ikb - 1 
    701                ze3 = fse3t_n(ji,jj,jk) 
    702                IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 
    703                IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) & 
    704                   &                                                       * r1_hisf_tbl(ji,jj) * ze3 
    705                IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) & 
    706                   &                                                       * r1_hisf_tbl(ji,jj) * ze3 
    707             END DO 
     705               DO jk = ikt, ikb - 1 
     706                  ze3 = fse3u_n(ji,jj,jk) 
     707                  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
     708               END DO 
    708709 
    709710            ! level partially include in ice shelf boundary layer  
    710             zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 
    711             IF (cptin == 'T') & 
    712                 &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 
    713             IF (cptin == 'U') & 
    714                 &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk) 
    715             IF (cptin == 'V') & 
    716                 &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk) 
    717          END DO 
    718       END DO 
     711               zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
     712               varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 
     713            END DO 
     714         END DO 
     715         DO jj = 2,jpj 
     716            DO ji = 2,jpi 
     717               varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji-1,jj)) 
     718            END DO 
     719         END DO 
     720         CALL lbc_lnk(varout,'T',-1.) 
     721      END IF 
     722 
     723      IF (cptin == 'V') THEN 
     724         DO jj = 1,jpj 
     725            DO ji = 1,jpi 
     726               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 
     727           ! thickness of boundary layer at least the top level thickness 
     728               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt)) 
     729 
     730            ! determine the deepest level influenced by the boundary layer 
     731            ! test on tmask useless ????? 
     732               DO jk = ikt, mbkv(ji,jj) 
     733                  IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 
     734               END DO 
     735               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     736 
     737            ! level fully include in the ice shelf boundary layer 
     738               DO jk = ikt, ikb - 1 
     739                  ze3 = fse3v_n(ji,jj,jk) 
     740                  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
     741               END DO 
     742 
     743            ! level partially include in ice shelf boundary layer  
     744               zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
     745               varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 
     746            END DO 
     747         END DO 
     748         DO jj = 2,jpj 
     749            DO ji = 2,jpi 
     750               varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji,jj-1)) 
     751            END DO 
     752         END DO 
     753         CALL lbc_lnk(varout,'T',-1.) 
     754      END IF 
     755 
     756      IF (cptin == 'T') THEN 
     757         DO jj = 1,jpj 
     758            DO ji = 1,jpi 
     759               ikt = misfkt(ji,jj) 
     760               ikb = misfkb(ji,jj) 
     761 
     762            ! level fully include in the ice shelf boundary layer 
     763               DO jk = ikt, ikb - 1 
     764                  ze3 = fse3t_n(ji,jj,jk) 
     765                  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 
     766               END DO 
     767 
     768            ! level partially include in ice shelf boundary layer  
     769               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 
     770               varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 
     771            END DO 
     772         END DO 
     773      END IF 
    719774      varout(:,:) = varout(:,:) * ssmask(:,:) 
    720775 
    721       CALL wrk_dealloc( jpi,jpj, mkt, mkb )       
    722       CALL wrk_dealloc( jpi,jpj, zikt )  
    723  
    724       IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.) 
    725       IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.) 
     776      CALL wrk_dealloc( jpi,jpj, zhisf_tbl )       
    726777 
    727778   END SUBROUTINE sbc_isf_tbl 
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r5204 r5216  
    271271               iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    272272               ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
    273                ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku) 
    274                ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv) 
     273               ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku) 
     274               ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv) 
    275275               ! 
    276276               ! i- direction 
     
    319319               iku = mbku(ji,jj) 
    320320               ikv = mbkv(ji,jj) 
    321                ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku) 
    322                ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv) 
    323  
    324                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1 
    325                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2 
    326                ENDIF 
    327                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1 
    328                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2 
     321               ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku) 
     322               ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv) 
     323 
     324               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)    ! i-direction: case 1 
     325               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)    ! -     -      case 2 
     326               ENDIF 
     327               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)    ! j-direction: case 1 
     328               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)    ! -     -      case 2 
    329329               ENDIF 
    330330 
     
    343343               iku = mbku(ji,jj) 
    344344               ikv = mbkv(ji,jj) 
    345                ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku) 
    346                ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv) 
     345               ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku) 
     346               ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv) 
    347347 
    348348               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1 
     
    370370               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 
    371371               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 
    372                ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku) 
    373                ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv)  
     372               ze3wu  =  fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku) 
     373               ze3wv  =  fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv)  
    374374 
    375375               ! i- direction 
     
    417417               iku = miku(ji,jj) 
    418418               ikv = mikv(ji,jj) 
    419                ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku) 
    420                ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv)  
    421  
    422                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1 
    423                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2 
    424                ENDIF 
    425  
    426                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1 
    427                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2 
     419               ze3wu  =  fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku) 
     420               ze3wv  =  fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv)  
     421 
     422               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)    ! i-direction: case 1 
     423               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)    ! -     -      case 2 
     424               ENDIF 
     425 
     426               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)    ! j-direction: case 1 
     427               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)    ! -     -      case 2 
    428428               ENDIF 
    429429 
     
    442442               iku = miku(ji,jj)  
    443443               ikv = mikv(ji,jj)  
    444                ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku) 
    445                ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv)  
     444               ze3wu  =  fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku) 
     445               ze3wv  =  fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv)  
    446446 
    447447               IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) ) ! i: 1 
Note: See TracChangeset for help on using the changeset viewer.