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

Changeset 5944


Ignore:
Timestamp:
2015-11-29T20:28:41+01:00 (8 years ago)
Author:
mathiot
Message:

ISF: change related to reviewers comments

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

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90

    r5621 r5944  
    137137      DO jk=1,nb_ana 
    138138       DO ji=1,jpmax_harmo 
    139           IF (TRIM(tname(jk)) .eq. Wave(ji)%cname_tide) THEN 
     139          IF (TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN 
    140140             name(jk) = ji 
    141141             EXIT 
     
    362362               X1=ana_amp(ji,jj,jh,1) 
    363363               X2=-ana_amp(ji,jj,jh,2) 
    364                out_v(ji,jj,       jh)=X1 * vmask_i(ji,jj) 
    365                out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj) 
     364               out_v(ji,jj,       jh)=X1 * ssvmask(ji,jj) 
     365               out_v(ji,jj,nb_ana+jh)=X2 * ssvmask(ji,jj) 
    366366            END DO 
    367367         END DO 
     
    492492            DO jj_sd = ji_sd, ninco 
    493493               zval2 = ABS(ztmp3(ji_sd,jj_sd)) 
    494                IF( zval2.GE.zval1 )THEN 
     494               IF( zval2 >= zval1 )THEN 
    495495                  ipivot(ji_sd) = jj_sd 
    496496                  zval1         = zval2 
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r5624 r5944  
    211211        CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content variation (psu) 
    212212        CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content variation (1.e20 J)  
    213         CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content variation (psu*km3) 
    214         CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh variation (km3)   
     213        CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9   )              ! Salt content variation (psu*km3) 
     214        CALL iom_put( 'bgvolssh' , zdiff_v1  * 1.e-9   )              ! volume ssh variation (km3)   
    215215        CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3)  
    216216        CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C)  
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r5621 r5944  
    319319         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla 
    320320      ENDIF 
    321       IF ( nn_cla .EQ. 1 ) THEN 
     321      IF ( nn_cla == 1 ) THEN 
    322322         IF  ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2  
    323323            CONTINUE 
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5624 r5944  
    332332         ! -------------------------------------------------- 
    333333         IF( ln_vvl_ztilde ) THEN 
    334             IF( kt .GT. nit000 ) THEN 
     334            IF( kt > nit000 ) THEN 
    335335               DO jk = 1, jpkm1 
    336336                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   & 
     
    426426         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
    427427         ! - ML - test: for the moment, stop simulation for too large e3_t variations 
    428          IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT. - rn_zdef_max ) ) THEN 
     428         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 
    429429            IF( lk_mpp ) THEN 
    430430               CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) 
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5621 r5944  
    824824   END SUBROUTINE zgr_bot_level 
    825825 
    826       SUBROUTINE zgr_top_level 
    827       !!---------------------------------------------------------------------- 
    828       !!                    ***  ROUTINE zgr_bot_level  *** 
     826   SUBROUTINE zgr_top_level 
     827      !!---------------------------------------------------------------------- 
     828      !!                    ***  ROUTINE zgr_top_level  *** 
    829829      !! 
    830830      !! ** Purpose :   defines the vertical index of ocean top (mik. arrays) 
     
    11491149                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    11501150               END DO 
    1151                IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 
     1151               IF (misfdep(ji,jj) >= 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 
    11521152               DO jk = misfdep(ji,jj) + 1, jpk 
    11531153                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     
    12121212      INTEGER  ::   ji, jj, jl, jk       ! dummy loop indices 
    12131213      INTEGER  ::   ik, it               ! temporary integers 
    1214       INTEGER  ::   icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF) 
     1214      INTEGER  ::   icompt, ibtest       ! (ISF) 
     1215      INTEGER  ::   ibtestim1, ibtestip1 ! (ISF) 
     1216      INTEGER  ::   ibtestjm1, ibtestjp1 ! (ISF) 
    12151217      REAL(wp) ::   zdepth           ! Ajusted ocean depth to avoid too small e3t 
    12161218      REAL(wp) ::   zmax             ! Maximum and minimum depth 
    1217       REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar 
     1219      REAL(wp) ::   zbathydiff       ! isf temporary scalar 
     1220      REAL(wp) ::   zrisfdepdiff     ! isf temporary scalar 
    12181221      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
    12191222      REAL(wp) ::   zdepwp           ! Ajusted ocean depth to avoid too small e3t 
     
    12301233 
    12311234      ! (ISF) compute misfdep 
    1232       WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
    1233       ELSEWHERE                      ;                          misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level  
     1235      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
     1236      ELSEWHERE                      ;                         misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level  
    12341237      END WHERE   
    12351238 
     
    12421245         WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1  
    12431246      END DO  
    1244       WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp) 
     1247      WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) ) 
    12451248         risfdep(:,:) = 0. ; misfdep(:,:) = 1 
    12461249      END WHERE 
    12471250 
    12481251      ! remove very shallow ice shelf (less than ~ 10m if 75L) 
    1249       IF ( cp_cfg .NE. "isomip" ) THEN 
    1250          WHERE (misfdep(:,:) <= 10 .AND. misfdep(:,:) .GT. 1) 
    1251             misfdep = 0; risfdep = 0.0_wp; 
    1252             mbathy  = 0; bathy   = 0.0_wp; 
    1253          END WHERE 
    1254          WHERE (bathy(:,:) <= 30.0 .AND. gphit < -60) 
    1255             misfdep = 0; risfdep = 0.0_wp; 
    1256             mbathy  = 0; bathy   = 0.0_wp; 
    1257          END WHERE 
    1258       END IF 
     1252      WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1) 
     1253         misfdep = 0; risfdep = 0.0_wp; 
     1254         mbathy  = 0; bathy   = 0.0_wp; 
     1255      END WHERE 
     1256      WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp) 
     1257         misfdep = 0; risfdep = 0.0_wp; 
     1258         mbathy  = 0; bathy   = 0.0_wp; 
     1259      END WHERE 
    12591260  
    12601261! 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 
     
    12621263! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 
    12631264      DO jl = 1, 10      
    1264          WHERE (bathy(:,:) .EQ. risfdep(:,:) ) 
     1265         WHERE (bathy(:,:) == risfdep(:,:) ) 
    12651266            misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
    12661267            mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
     
    12711272         ENDWHERE 
    12721273         IF( lk_mpp ) THEN 
    1273             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1274            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    12741275            CALL lbc_lnk( zbathy, 'T', 1. ) 
    12751276            misfdep(:,:) = INT( zbathy(:,:) ) 
    1276             CALL lbc_lnk( risfdep, 'T', 1. ) 
    1277             CALL lbc_lnk( bathy, 'T', 1. ) 
    1278             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1277 
     1278            CALL lbc_lnk( risfdep,'T', 1. ) 
     1279            CALL lbc_lnk( bathy,  'T', 1. ) 
     1280 
     1281            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    12791282            CALL lbc_lnk( zbathy, 'T', 1. ) 
    1280             mbathy(:,:) = INT( zbathy(:,:) ) 
     1283            mbathy(:,:)  = INT( zbathy(:,:) ) 
    12811284         ENDIF 
    12821285         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN  
    1283             misfdep( 1 ,:) = misfdep(jpim1,:)           ! local domain is cyclic east-west  
     1286            misfdep( 1 ,:) = misfdep(jpim1,:)            ! local domain is cyclic east-west  
    12841287            misfdep(jpi,:) = misfdep(  2  ,:)  
    1285          ENDIF 
    1286  
    1287          IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
    1288             mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west 
    1289             mbathy(jpi,:) = mbathy(  2  ,:) 
     1288            mbathy( 1 ,:)  = mbathy(jpim1,:)             ! local domain is cyclic east-west 
     1289            mbathy(jpi,:)  = mbathy(  2  ,:) 
    12901290         ENDIF 
    12911291 
     
    13141314               ! find the minimum change option: 
    13151315               ! test bathy 
    1316                IF (risfdep(ji,jj) .GT. 1) THEN 
     1316               IF (risfdep(ji,jj) > 1) THEN 
    13171317               zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) & 
    13181318                 &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
     
    13201320                 &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    13211321  
    1322                   IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 
    1323                      IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1322                  IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 
     1323                     IF (zbathydiff <= zrisfdepdiff) THEN 
    13241324                        bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
    13251325                        mbathy(ji,jj)= mbathy(ji,jj) + 1 
     
    13381338               ! find the minimum change option: 
    13391339               ! test bathy 
    1340                IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1340               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 
    13411341                  zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1)& 
    13421342                    & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
    13431343                  zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  )  & 
    13441344                    & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    1345                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1345                  IF (zbathydiff <= zrisfdepdiff) THEN 
    13461346                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
    13471347                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
     
    13541354         END DO 
    13551355  
    1356  ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)  
     1356 ! point V mbathy(ji,jj) == misfdep(ji,jj+1)  
    13571357         DO jj = 1, jpjm1 
    13581358            DO ji = 1, jpim1 
    1359                IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1359               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 
    13601360                  zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1)   & 
    13611361                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat ))) 
    13621362                  zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1))  & 
    13631363                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 
    1364                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1364                  IF (zbathydiff <= zrisfdepdiff) THEN 
    13651365                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
    13661366                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) & 
     
    13761376  
    13771377         IF( lk_mpp ) THEN 
    1378             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1378            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    13791379            CALL lbc_lnk( zbathy, 'T', 1. ) 
    13801380            misfdep(:,:) = INT( zbathy(:,:) ) 
    1381             CALL lbc_lnk( risfdep, 'T', 1. ) 
    1382             CALL lbc_lnk( bathy, 'T', 1. ) 
    1383             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1381 
     1382            CALL lbc_lnk( risfdep,'T', 1. ) 
     1383            CALL lbc_lnk( bathy,  'T', 1. ) 
     1384 
     1385            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    13841386            CALL lbc_lnk( zbathy, 'T', 1. ) 
    1385             mbathy(:,:) = INT( zbathy(:,:) ) 
     1387            mbathy(:,:)  = INT( zbathy(:,:) ) 
    13861388         ENDIF 
    1387  ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)  
     1389 ! point V misdep(ji,jj) == mbathy(ji,jj+1)  
    13881390         DO jj = 1, jpjm1 
    13891391            DO ji = 1, jpim1 
    1390                IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1392               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN 
    13911393                  zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) & 
    13921394                   &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 
    13931395                  zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  )  & 
    13941396                   &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat ))) 
    1395                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1397                  IF (zbathydiff <= zrisfdepdiff) THEN 
    13961398                     mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 
    13971399                     bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 
     
    14061408  
    14071409         IF( lk_mpp ) THEN  
    1408             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1409             CALL lbc_lnk( zbathy, 'T', 1. )  
    1410             misfdep(:,:) = INT( zbathy(:,:) )  
    1411             CALL lbc_lnk( risfdep, 'T', 1. )  
    1412             CALL lbc_lnk( bathy, 'T', 1. ) 
    1413             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1410            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    14141411            CALL lbc_lnk( zbathy, 'T', 1. ) 
    1415             mbathy(:,:) = INT( zbathy(:,:) ) 
     1412            misfdep(:,:) = INT( zbathy(:,:) ) 
     1413 
     1414            CALL lbc_lnk( risfdep,'T', 1. ) 
     1415            CALL lbc_lnk( bathy,  'T', 1. ) 
     1416 
     1417            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
     1418            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1419            mbathy(:,:)  = INT( zbathy(:,:) ) 
    14161420         ENDIF  
    14171421  
    1418  ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)  
     1422 ! point U mbathy(ji,jj) == misfdep(ji,jj+1)  
    14191423         DO jj = 1, jpjm1 
    14201424            DO ji = 1, jpim1 
    1421                IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1425               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 
    14221426                  zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 
    14231427                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat ))) 
    14241428                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) & 
    14251429                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 
    1426                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1430                  IF (zbathydiff <= zrisfdepdiff) THEN 
    14271431                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
    14281432                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
     
    14361440  
    14371441         IF( lk_mpp ) THEN  
    1438             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1439             CALL lbc_lnk( zbathy, 'T', 1. )  
    1440             misfdep(:,:) = INT( zbathy(:,:) )  
    1441             CALL lbc_lnk( risfdep, 'T', 1. )  
    1442             CALL lbc_lnk( bathy, 'T', 1. ) 
    1443             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1442            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    14441443            CALL lbc_lnk( zbathy, 'T', 1. ) 
    1445             mbathy(:,:) = INT( zbathy(:,:) ) 
     1444            misfdep(:,:) = INT( zbathy(:,:) ) 
     1445 
     1446            CALL lbc_lnk( risfdep,'T', 1. ) 
     1447            CALL lbc_lnk( bathy,  'T', 1. ) 
     1448 
     1449            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
     1450            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1451            mbathy(:,:)  = INT( zbathy(:,:) ) 
    14461452         ENDIF  
    14471453  
    1448  ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)  
     1454 ! point U misfdep(ji,jj) == bathy(ji,jj+1)  
    14491455         DO jj = 1, jpjm1 
    14501456            DO ji = 1, jpim1 
    1451                IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
     1457               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN 
    14521458                  zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) & 
    14531459                      &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 
    14541460                  zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  )  & 
    14551461                      &  - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat ))) 
    1456                   IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1462                  IF (zbathydiff <= zrisfdepdiff) THEN 
    14571463                     mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1 
    14581464                     bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  )  & 
     
    14681474  
    14691475         IF( lk_mpp ) THEN 
    1470             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1476            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    14711477            CALL lbc_lnk( zbathy, 'T', 1. ) 
    14721478            misfdep(:,:) = INT( zbathy(:,:) ) 
    1473             CALL lbc_lnk( risfdep, 'T', 1. ) 
    1474             CALL lbc_lnk( bathy, 'T', 1. ) 
    1475             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1479 
     1480            CALL lbc_lnk( risfdep,'T', 1. ) 
     1481            CALL lbc_lnk( bathy,  'T', 1. ) 
     1482 
     1483            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    14761484            CALL lbc_lnk( zbathy, 'T', 1. ) 
    1477             mbathy(:,:) = INT( zbathy(:,:) ) 
     1485            mbathy(:,:)  = INT( zbathy(:,:) ) 
    14781486         ENDIF 
    14791487      END DO 
     
    14851493         DO jk = 2, jpk 
    14861494            WHERE (misfdep==0) misfdep=jpk 
    1487             zmask=0 
    1488             WHERE (misfdep .LE. jk) zmask=1 
     1495            zmask=0._wp 
     1496            WHERE (misfdep <= jk) zmask=1 
    14891497            DO jj = 2, jpjm1 
    14901498               DO ji = 2, jpim1 
    1491                   IF (misfdep(ji,jj) .EQ. jk) THEN 
     1499                  IF (misfdep(ji,jj) == jk) THEN 
    14921500                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
    1493                      IF (ibtest .LE. 1) THEN 
     1501                     IF (ibtest <= 1) THEN 
    14941502                        risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 
    1495                         IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk 
     1503                        IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk 
    14961504                     END IF 
    14971505                  END IF 
     
    15001508         END DO 
    15011509         WHERE (misfdep==jpk) 
    1502              misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 
     1510             misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 
    15031511         END WHERE 
    15041512         IF( lk_mpp ) THEN 
    1505             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1513            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    15061514            CALL lbc_lnk( zbathy, 'T', 1. ) 
    15071515            misfdep(:,:) = INT( zbathy(:,:) ) 
    1508             CALL lbc_lnk( risfdep, 'T', 1. ) 
    1509             CALL lbc_lnk( bathy, 'T', 1. ) 
    1510             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1516 
     1517            CALL lbc_lnk( risfdep,'T', 1. ) 
     1518            CALL lbc_lnk( bathy,  'T', 1. ) 
     1519 
     1520            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    15111521            CALL lbc_lnk( zbathy, 'T', 1. ) 
    1512             mbathy(:,:) = INT( zbathy(:,:) ) 
     1522            mbathy(:,:)  = INT( zbathy(:,:) ) 
    15131523         ENDIF 
    15141524  
    15151525 ! remove single point "bay" on bathy coast line beneath an ice shelf' 
    15161526         DO jk = jpk,1,-1 
    1517             zmask=0 
    1518             WHERE (mbathy .GE. jk ) zmask=1 
     1527            zmask=0._wp 
     1528            WHERE (mbathy >= jk ) zmask=1 
    15191529            DO jj = 2, jpjm1 
    15201530               DO ji = 2, jpim1 
    1521                   IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN 
     1531                  IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN 
    15221532                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
    1523                      IF (ibtest .LE. 1) THEN 
     1533                     IF (ibtest <= 1) THEN 
    15241534                        bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 
    1525                         IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0 
     1535                        IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0 
    15261536                     END IF 
    15271537                  END IF 
     
    15301540         END DO 
    15311541         WHERE (mbathy==0) 
    1532              misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 
     1542             misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 
    15331543         END WHERE 
    15341544         IF( lk_mpp ) THEN 
    1535             zbathy(:,:) = FLOAT( misfdep(:,:) ) 
     1545            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    15361546            CALL lbc_lnk( zbathy, 'T', 1. ) 
    15371547            misfdep(:,:) = INT( zbathy(:,:) ) 
    1538             CALL lbc_lnk( risfdep, 'T', 1. ) 
    1539             CALL lbc_lnk( bathy, 'T', 1. ) 
    1540             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1548 
     1549            CALL lbc_lnk( risfdep,'T', 1. ) 
     1550            CALL lbc_lnk( bathy,  'T', 1. ) 
     1551 
     1552            zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    15411553            CALL lbc_lnk( zbathy, 'T', 1. ) 
    1542             mbathy(:,:) = INT( zbathy(:,:) ) 
     1554            mbathy(:,:)  = INT( zbathy(:,:) ) 
    15431555         ENDIF 
    15441556  
     
    15461558         zmisfdep = misfdep 
    15471559         zrisfdep = risfdep 
    1548          WHERE (zmisfdep .LE. 1) zmisfdep=jpk 
     1560         WHERE (zmisfdep <= 1._wp) zmisfdep=jpk 
    15491561         DO jj = 2, jpjm1 
    15501562            DO ji = 2, jpim1 
    15511563               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  ) 
    15521564               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1) 
    1553                IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1) 
    1554                IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1) 
    1555                IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1) 
    1556                IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1) 
     1565               IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1) 
     1566               IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1) 
     1567               IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1) 
     1568               IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1) 
    15571569               ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1558                IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN 
     1570               IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN 
    15591571                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 
    15601572               END IF 
    1561                IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN 
     1573               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN 
    15621574                  misfdep(ji,jj) = ibtest 
    15631575                  risfdep(ji,jj) = gdepw_1d(ibtest) 
     
    15671579  
    15681580         IF( lk_mpp ) THEN  
    1569             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1570             CALL lbc_lnk( zbathy, 'T', 1. )  
     1581            zbathy(:,:)  = FLOAT( misfdep(:,:) )  
     1582            CALL lbc_lnk( zbathy,  'T', 1. )  
    15711583            misfdep(:,:) = INT( zbathy(:,:) )  
     1584 
    15721585            CALL lbc_lnk( risfdep, 'T', 1. )  
    1573             CALL lbc_lnk( bathy, 'T', 1. ) 
     1586            CALL lbc_lnk( bathy,   'T', 1. ) 
     1587 
    15741588            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1575             CALL lbc_lnk( zbathy, 'T', 1. ) 
     1589            CALL lbc_lnk( zbathy,  'T', 1. ) 
    15761590            mbathy(:,:) = INT( zbathy(:,:) ) 
    15771591         ENDIF  
     
    15831597               ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  ) 
    15841598               ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1) 
    1585                IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1) 
    1586                IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj  ) ) ibtestip1 = 0 
    1587                IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj-1) ) ibtestjm1 = 0 
    1588                IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
     1599               IF( zmbathy(ji,jj) < misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1) 
     1600               IF( zmbathy(ji,jj) < misfdep(ji+1,jj  ) ) ibtestip1 = 0 
     1601               IF( zmbathy(ji,jj) < misfdep(ji  ,jj-1) ) ibtestjm1 = 0 
     1602               IF( zmbathy(ji,jj) < misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
    15891603               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1590                IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 
     1604               IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN 
    15911605                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
    15921606               END IF 
    1593                IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN 
     1607               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN 
    15941608                  mbathy(ji,jj) = ibtest 
    15951609                  bathy(ji,jj)  = gdepw_1d(ibtest+1)  
     
    15981612         END DO 
    15991613         IF( lk_mpp ) THEN  
    1600             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1601             CALL lbc_lnk( zbathy, 'T', 1. )  
     1614            zbathy(:,:)  = FLOAT( misfdep(:,:) )  
     1615            CALL lbc_lnk( zbathy,  'T', 1. )  
    16021616            misfdep(:,:) = INT( zbathy(:,:) )  
     1617 
    16031618            CALL lbc_lnk( risfdep, 'T', 1. )  
    1604             CALL lbc_lnk( bathy, 'T', 1. ) 
     1619            CALL lbc_lnk( bathy,   'T', 1. ) 
     1620 
    16051621            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1606             CALL lbc_lnk( zbathy, 'T', 1. ) 
     1622            CALL lbc_lnk( zbathy,  'T', 1. ) 
    16071623            mbathy(:,:) = INT( zbathy(:,:) ) 
    16081624         ENDIF  
     
    16101626         DO jj = 1, jpjm1 
    16111627            DO ji = 1, jpim1 
    1612                IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 
     1628               IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 
    16131629                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
    16141630               END IF 
     
    16161632         END DO 
    16171633         IF( lk_mpp ) THEN  
    1618             zbathy(:,:) = FLOAT( misfdep(:,:) )  
    1619             CALL lbc_lnk( zbathy, 'T', 1. )  
     1634            zbathy(:,:)  = FLOAT( misfdep(:,:) )  
     1635            CALL lbc_lnk( zbathy,  'T', 1. )  
    16201636            misfdep(:,:) = INT( zbathy(:,:) )  
     1637 
    16211638            CALL lbc_lnk( risfdep, 'T', 1. )  
    1622             CALL lbc_lnk( bathy, 'T', 1. ) 
     1639            CALL lbc_lnk( bathy,   'T', 1. ) 
     1640 
    16231641            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1624             CALL lbc_lnk( zbathy, 'T', 1. ) 
     1642            CALL lbc_lnk( zbathy,  'T', 1. ) 
    16251643            mbathy(:,:) = INT( zbathy(:,:) ) 
    16261644         ENDIF  
     
    16281646         DO jj = 1, jpjm1 
    16291647            DO ji = 1, jpim1 
    1630                IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 
     1648               IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 
    16311649                  mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ; 
    16321650               END IF 
     
    16341652         END DO 
    16351653         IF( lk_mpp ) THEN  
    1636             zbathy(:,:) = FLOAT( misfdep(:,:) )  
     1654            zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    16371655            CALL lbc_lnk( zbathy, 'T', 1. )  
    16381656            misfdep(:,:) = INT( zbathy(:,:) )  
    1639             CALL lbc_lnk( risfdep, 'T', 1. )  
    1640             CALL lbc_lnk( bathy, 'T', 1. ) 
     1657 
     1658            CALL lbc_lnk( risfdep,'T', 1. )  
     1659            CALL lbc_lnk( bathy,  'T', 1. ) 
     1660 
    16411661            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    16421662            CALL lbc_lnk( zbathy, 'T', 1. ) 
     
    16461666         DO jj = 1, jpjm1 
    16471667            DO ji = 1, jpi 
    1648                IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 
     1668               IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 
    16491669                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
    16501670               END IF 
     
    16521672         END DO 
    16531673         IF( lk_mpp ) THEN  
    1654             zbathy(:,:) = FLOAT( misfdep(:,:) )  
     1674            zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    16551675            CALL lbc_lnk( zbathy, 'T', 1. )  
    16561676            misfdep(:,:) = INT( zbathy(:,:) )  
    1657             CALL lbc_lnk( risfdep, 'T', 1. )  
    1658             CALL lbc_lnk( bathy, 'T', 1. ) 
     1677 
     1678            CALL lbc_lnk( risfdep,'T', 1. )  
     1679            CALL lbc_lnk( bathy,  'T', 1. ) 
     1680 
    16591681            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    16601682            CALL lbc_lnk( zbathy, 'T', 1. ) 
     
    16641686         DO jj = 1, jpjm1 
    16651687            DO ji = 1, jpi 
    1666                IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 
    1667                   mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1)   = gdepw_1d(mbathy(ji,jj+1)+1) ; 
     1688               IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 
     1689                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 
    16681690               END IF 
    16691691            END DO 
    16701692         END DO 
    16711693         IF( lk_mpp ) THEN  
    1672             zbathy(:,:) = FLOAT( misfdep(:,:) )  
     1694            zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    16731695            CALL lbc_lnk( zbathy, 'T', 1. )  
    16741696            misfdep(:,:) = INT( zbathy(:,:) )  
    1675             CALL lbc_lnk( risfdep, 'T', 1. )  
    1676             CALL lbc_lnk( bathy, 'T', 1. ) 
     1697 
     1698            CALL lbc_lnk( risfdep,'T', 1. )  
     1699            CALL lbc_lnk( bathy,  'T', 1. ) 
     1700 
    16771701            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    16781702            CALL lbc_lnk( zbathy, 'T', 1. ) 
     
    16941718! end check compatibility ice shelf/bathy 
    16951719      ! remove very shallow ice shelf (less than ~ 10m if 75L) 
    1696       WHERE (misfdep(:,:) <= 5) 
     1720      WHERE (risfdep(:,:) <= 10._wp) 
    16971721         misfdep = 1; risfdep = 0.0_wp; 
    16981722      END WHERE 
     
    18141838      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf') 
    18151839       
    1816    END SUBROUTINE 
     1840   END SUBROUTINE zgr_isf 
    18171841 
    18181842   SUBROUTINE zgr_sco 
     
    21432167      fsde3w(:,:,:) = gdep3w_0(:,:,:) 
    21442168      ! 
    2145       where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0 
    2146       where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0 
    2147       where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0 
    2148       where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0 
    2149       where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0 
    2150       where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0 
    2151       where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0 
     2169      where (e3t_0   (:,:,:) == 0.0)  e3t_0(:,:,:)  = 1.0_wp 
     2170      where (e3u_0   (:,:,:) == 0.0)  e3u_0(:,:,:)  = 1.0_wp 
     2171      where (e3v_0   (:,:,:) == 0.0)  e3v_0(:,:,:)  = 1.0_wp 
     2172      where (e3f_0   (:,:,:) == 0.0)  e3f_0(:,:,:)  = 1.0_wp 
     2173      where (e3w_0   (:,:,:) == 0.0)  e3w_0(:,:,:)  = 1.0_wp 
     2174      where (e3uw_0  (:,:,:) == 0.0)  e3uw_0(:,:,:) = 1.0_wp 
     2175      where (e3vw_0  (:,:,:) == 0.0)  e3vw_0(:,:,:) = 1.0_wp 
    21522176 
    21532177#if defined key_agrif 
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r5921 r5944  
    235235                     &                              * (1._wp - tmask(ji,jj,jk)) 
    236236               END DO 
    237                IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 
     237               IF (ikt  >= 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 
    238238                                  &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 
    239239            END DO 
     
    521521      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj 
    522522      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop 
    523       REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_oce, zpshpi, zpshpj 
     523      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_oce 
    524524      !!---------------------------------------------------------------------- 
    525525      ! 
    526526      CALL wrk_alloc( jpi,jpj,  2, ztstop)  
    527527      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj) 
    528       CALL wrk_alloc( jpi,jpj,     zrhdtop_oce, zpshpi, zpshpj)  
     528      CALL wrk_alloc( jpi,jpj,     zrhdtop_oce ) 
    529529      ! 
    530530      ! Local constant initialization 
     
    605605      CALL wrk_dealloc( jpi,jpj,2  , ztstop) 
    606606      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj) 
    607       CALL wrk_dealloc( jpi,jpj    , zrhdtop_oce, zpshpi, zpshpj) 
     607      CALL wrk_dealloc( jpi,jpj    , zrhdtop_oce ) 
    608608      ! 
    609609   END SUBROUTINE hpg_isf 
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r5624 r5944  
    100100      ! 
    101101      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    102       INTEGER  ::   iku, ikv     ! local integers 
     102      INTEGER  ::   ikt          ! local integers 
    103103#if ! defined key_dynspg_flt 
    104104      REAL(wp) ::   z2dt         ! temporary scalar 
     
    268268               DO jj = 1,jpj 
    269269                  DO ji = 1,jpi 
    270                      jk = mikt(ji,jj) 
    271                      fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0                  &  
     270                     ikt = mikt(ji,jj) 
     271                     fse3t_b(ji,jj,ikt) = fse3t_b(ji,jj,ikt) - atfp * rdt * r1_rau0                  &  
    272272                                       &                          * ( (emp_b(ji,jj) - emp(ji,jj) ) & 
    273273                                       &                            - (rnf_b(ji,jj) - rnf(ji,jj) ) &     
    274                                        &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj) ) ) * tmask(ji,jj,jk) 
     274                                       &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj) ) ) * tmask(ji,jj,ikt) 
    275275                  END DO 
    276276               END DO 
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r5910 r5944  
    934934      ! 
    935935      IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN 
    936          IF ( Agrif_NbStepint().EQ.0 ) THEN 
     936         IF ( Agrif_NbStepint() == 0 ) THEN 
    937937            ub2_i_b(:,:) = 0.e0 
    938938            vb2_i_b(:,:) = 0.e0 
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

    r5189 r5944  
    8686         DO jj = 2, jpj                   ! vertical fluxes  
    8787            DO ji = fs_2, jpi 
    88                zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
     88               zww(ji,jj) = 0.25_wp * e12t(ji,jj) * wn(ji,jj,jk) 
    8989            END DO 
    9090         END DO 
     
    193193         DO jj = 2, jpj                    
    194194            DO ji = fs_2, jpi             ! vector opt. 
    195                zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 
     195               zww(ji,jj,jk) = 0.25_wp * e12t(ji,jj) * wn(ji,jj,jk) 
    196196            END DO 
    197197         END DO 
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r5200 r5944  
    586586                  ! 
    587587                  jk = nmln(ji+ip,jj) + 1 
    588                   IF( jk .GT. mbkt(ji+ip,jj) ) THEN  !ML reaches bottom 
     588                  IF( jk > mbkt(ji+ip,jj) ) THEN  !ML reaches bottom 
    589589                    zti_mlb(ji+ip,jj   ,1-ip,kp) = 0.0_wp 
    590590                  ELSE 
     
    596596                  ! 
    597597                  jk = nmln(ji,jj+jp) + 1 
    598                   IF( jk .GT. mbkt(ji,jj+jp) ) THEN  !ML reaches bottom 
     598                  IF( jk > mbkt(ji,jj+jp) ) THEN  !ML reaches bottom 
    599599                    ztj_mlb(ji   ,jj+jp,1-jp,kp) = 0.0_wp 
    600600                  ELSE 
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90

    r5624 r5944  
    9191         IF( kn_fwb == 3 .AND. ln_isfcav    )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 with ln_isfcav = .TRUE. not working, we stop ' ) 
    9292         ! 
    93          area = glob_sum( e1e2t(:,:) * tmask(:,:,1))           ! interior global domain surface 
     93         area = glob_sum( e12t(:,:) * tmask(:,:,1))           ! interior global domain surface 
    9494         ! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes 
    9595         ! and in case of no melt, it can generate HSSW. 
     
    108108         ! 
    109109         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN 
    110             z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain 
     110            z_fwf = glob_sum( e12t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain 
    111111            zcoef = z_fwf * rcp 
    112112            emp(:,:) = emp(:,:) - z_fwf              * tmask(:,:,1) 
     
    133133            a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow 
    134134                                                      ! sum over the global domain 
    135             a_fwb   = glob_sum( e1e2t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) ) 
     135            a_fwb   = glob_sum( e12t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) ) 
    136136            a_fwb   = a_fwb * 1.e+3 / ( area * rday * 365. )     ! convert in Kg/m3/s = mm/s 
    137137!!gm        !                                                      !!bug 365d year  
     
    159159            ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:) 
    160160            ! 
    161             zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )  ! Area filled by <0 and >0 erp  
    162             zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) ) 
     161            zsurf_neg = glob_sum( e12t(:,:)*ztmsk_neg(:,:) )  ! Area filled by <0 and >0 erp  
     162            zsurf_pos = glob_sum( e12t(:,:)*ztmsk_pos(:,:) ) 
    163163            !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges)  
    164             z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area 
     164            z_fwf     = glob_sum( e12t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area 
    165165            !             
    166166            IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation 
     
    172172            ENDIF 
    173173            ! 
    174             zsum_fwf   = glob_sum( e1e2t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area 
     174            zsum_fwf   = glob_sum( e12t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area 
    175175!!gm :  zsum_fwf   = z_fwf * area   ???  it is right?  I think so.... 
    176176            z_fwf_nsrf =  zsum_fwf / ( zsurf_tospread + rsmall ) 
    177177            !                                                  ! weight to respect erp field 2D structure  
    178             zsum_erp   = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) ) 
     178            zsum_erp   = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e12t(:,:) ) 
    179179            z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall ) 
    180180            !                                                  ! final correction term to apply 
     
    191191               IF( z_fwf < 0._wp ) THEN 
    192192                  WRITE(numout,*)'   z_fwf < 0' 
    193                   WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv' 
     193                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e12t(:,:) )*1.e-9,' Sv' 
    194194               ELSE 
    195195                  WRITE(numout,*)'   z_fwf >= 0' 
    196                   WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv' 
     196                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e12t(:,:) )*1.e-9,' Sv' 
    197197               ENDIF 
    198                WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e1e2t(:,:) )*1.e-9,' Sv' 
     198               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e12t(:,:) )*1.e-9,' Sv' 
    199199               WRITE(numout,*)'   z_fwf         = ', z_fwf      ,' Kg/m2/s' 
    200200               WRITE(numout,*)'   z_fwf_nsrf    = ', z_fwf_nsrf ,' Kg/m2/s' 
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r5905 r5944  
    5656                                                                                          !: (first wet level and last level include in the tbl) 
    5757#else 
    58    INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base 
     58   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)      ::  misfkt, misfkb         !:Level of ice shelf base 
    5959#endif 
    6060 
    6161 
    62    REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ? 
    63    REAL(wp), PUBLIC, SAVE ::   kappa  = 1.54e-6_wp    ! phycst ? 
    64    REAL(wp), PUBLIC, SAVE ::   rhoisf = 920.0_wp      ! phycst ? 
    65    REAL(wp), PUBLIC, SAVE ::   tsurf  = -20.0_wp      ! phycst ? 
    66    REAL(wp), PUBLIC, SAVE ::   lfusisf= 0.334e6_wp    ! phycst ? 
     62   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     ! phycst ? 
     63   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    ! phycst ? 
     64   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      ! phycst ? 
     65   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      ! phycst ? 
     66   REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    ! phycst ? 
    6767 
    6868!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 
    69    CHARACTER(len=100), PUBLIC ::   cn_dirisf  = './'    !: Root directory for location of ssr files 
    70    TYPE(FLD_N)       , PUBLIC ::   sn_fwfisf            !: information about the isf file to be read 
    71    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_fwfisf 
    72    TYPE(FLD_N)       , PUBLIC ::   sn_rnfisf              !: information about the runoff file to be read 
    73    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnfisf            
    74    TYPE(FLD_N)       , PUBLIC ::   sn_depmax_isf, sn_depmin_isf, sn_Leff_isf     !: information about the runoff file to be read 
     69   CHARACTER(len=100), PUBLIC           :: cn_dirisf  = './' !: Root directory for location of ssr files 
     70   TYPE(FLD_N)       , PUBLIC           :: sn_fwfisf         !: information about the isf file to be read 
     71   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf 
     72   TYPE(FLD_N)       , PUBLIC           :: sn_rnfisf         !: information about the runoff file to be read 
     73   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf            
     74   TYPE(FLD_N)       , PUBLIC           :: sn_depmax_isf     !: information about the runoff file to be read ?? 
     75   TYPE(FLD_N)       , PUBLIC           :: sn_depmin_isf     !: information about the runoff file to be read ?? 
     76   TYPE(FLD_N)       , PUBLIC           :: sn_Leff_isf       !: information about the runoff file to be read ?? 
    7577    
    7678   !! * Substitutions 
     
    8587  
    8688  SUBROUTINE sbc_isf(kt) 
    87     INTEGER, INTENT(in) :: kt                   ! ocean time step 
    88     INTEGER             :: ji, jj, jk           ! loop index 
    89     INTEGER             :: ikt, ikb             ! top and bottom level of the isf boundary layer 
    90     INTEGER             :: inum, ierror 
    91     INTEGER             :: ios                  ! Local integer output status for namelist read 
    92     REAL(wp)            :: zhk 
    93     REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep)  
    94     CHARACTER(len=256)  :: cvarzisf, cvarhisf   ! name for isf file 
    95     CHARACTER(LEN=32 )  :: cvarLeff             ! variable name for efficient Length scale 
    96       ! 
    9789      !!--------------------------------------------------------------------- 
    98       NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf      , & 
    99                          & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 
    100       ! 
    101       ! allocation 
    102       CALL wrk_alloc( jpi,jpj, zt_frz, zdep  ) 
     90      !!                  ***  ROUTINE sbc_isf  *** 
     91      !! 
     92      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf  
     93      !!              melting and freezing  
     94      !! 
     95      !! ** Method  :  4 parameterizations are available according to nn_isf  
     96      !!               nn_isf = 1 : Realistic ice_shelf formulation 
     97      !!                        2 : Beckmann & Goose parameterization 
     98      !!                        3 : Specified runoff in deptht (Mathiot & al. ) 
     99      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf 
     100      !!---------------------------------------------------------------------- 
     101      INTEGER, INTENT( in ) :: kt                   ! ocean time step 
     102      ! 
     103      INTEGER               :: ji, jj               ! loop index 
     104      REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep)  
     105      !!--------------------------------------------------------------------- 
    103106      ! 
    104107      !                                         ! ====================== ! 
    105108      IF( kt == nit000 ) THEN                   !  First call kt=nit000  ! 
    106109         !                                      ! ====================== ! 
    107          REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs  
    108          READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) 
    109 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp ) 
    110  
    111          REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs 
    112          READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 ) 
    113 902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 
    114          IF(lwm) WRITE ( numond, namsbc_isf ) 
    115  
    116  
    117          IF ( lwp ) WRITE(numout,*) 
    118          IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf' 
    119          IF ( lwp ) WRITE(numout,*) '~~~~~~~~~' 
    120          IF ( lwp ) WRITE(numout,*) 'sbcisf :'  
    121          IF ( lwp ) WRITE(numout,*) '~~~~~~~~' 
    122          IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf 
    123          IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk 
    124          IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl 
    125          IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk  
    126          IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0   
    127          IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0   
    128          IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2  
     110         CALL sbc_isf_init 
     111      !                                         ! ---------------------------------------- ! 
     112      ELSE                                      !          Swap of forcing fields          ! 
     113         !                                      ! ---------------------------------------- ! 
     114         fwfisf_b  (:,:  ) = fwfisf  (:,:  )    ! Swap the ocean forcing fields except at nit000 
     115         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)    ! where before fields are set at the end of the routine 
    129116         ! 
    130          ! Allocate public variable 
    131          IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 
    132          ! 
    133          ! initialisation 
    134          qisf(:,:)        = 0._wp  ; fwfisf  (:,:) = 0._wp 
    135          risf_tsc(:,:,:)  = 0._wp  ; fwfisf_b(:,:) = 0._wp 
    136          ! 
    137          ! define isf tbl tickness, top and bottom indice 
    138          IF      (nn_isf == 1) THEN 
    139             rhisf_tbl(:,:) = rn_hisf_tbl 
    140             misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv 
    141          ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN 
    142             ALLOCATE( sf_rnfisf(1), STAT=ierror ) 
    143             ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 
    144             CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
    145  
    146             !: read effective lenght (BG03) 
    147             IF (nn_isf == 2) THEN 
    148                cvarLeff = 'soLeff'  
    149                CALL iom_open( sn_Leff_isf%clname, inum ) 
    150                CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 
    151                CALL iom_close(inum) 
    152                ! 
    153                risfLeff = risfLeff*1000.0_wp           !: convertion in m 
    154             END IF 
    155  
    156            ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth) 
    157             CALL iom_open( sn_depmax_isf%clname, inum ) 
    158             cvarhisf = TRIM(sn_depmax_isf%clvar) 
    159             CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base 
    160             CALL iom_close(inum) 
    161             ! 
    162             CALL iom_open( sn_depmin_isf%clname, inum ) 
    163             cvarzisf = TRIM(sn_depmin_isf%clvar) 
    164             CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base 
    165             CALL iom_close(inum) 
    166             ! 
    167             rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer 
    168  
    169            !! compute first level of the top boundary layer 
    170            DO ji = 1, jpi 
    171               DO jj = 1, jpj 
    172                   jk = 2 
    173                   DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdepw(ji,jj,jk) < rzisf_tbl(ji,jj) ) ;  jk = jk + 1 ;  END DO 
    174                   misfkt(ji,jj) = jk-1 
    175                END DO 
    176             END DO 
    177  
    178          ELSE IF ( nn_isf == 4 ) THEN 
    179             ! as in nn_isf == 1 
    180             rhisf_tbl(:,:) = rn_hisf_tbl 
    181             misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv 
    182              
    183             ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 
    184             ALLOCATE( sf_fwfisf(1), STAT=ierror ) 
    185             ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 
    186             CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
    187          END IF 
    188           
    189          rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 
    190  
    191          ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 
    192          DO jj = 1,jpj 
    193             DO ji = 1,jpi 
    194                ikt = misfkt(ji,jj) 
    195                ikb = misfkt(ji,jj) 
    196                ! thickness of boundary layer at least the top level thickness 
    197                rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt)) 
    198  
    199                ! determine the deepest level influenced by the boundary layer 
    200                DO jk = ikt+1, mbkt(ji,jj) 
    201                   IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
    202                END DO 
    203                rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
    204                misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl 
    205                r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 
    206  
    207                zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1 
    208                ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer 
    209             END DO 
    210          END DO 
    211           
    212117      END IF 
    213118 
    214       !                                            ! ---------------------------------------- ! 
    215       IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          ! 
    216          !                                         ! ---------------------------------------- ! 
    217          fwfisf_b  (:,:  ) = fwfisf  (:,:  )               ! Swap the ocean forcing fields except at nit000 
    218          risf_tsc_b(:,:,:) = risf_tsc(:,:,:)               ! where before fields are set at the end of the routine 
    219          ! 
    220       ENDIF 
    221  
    222119      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN 
    223  
    224  
    225          ! compute salf and heat flux 
    226          IF (nn_isf == 1) THEN 
    227             ! realistic ice shelf formulation 
     120         ! allocation 
     121         CALL wrk_alloc( jpi,jpj, zt_frz, zdep  ) 
     122 
     123         ! compute salt and heat flux 
     124         SELECT CASE ( nn_isf ) 
     125         CASE ( 1 )    ! realistic ice shelf formulation 
    228126            ! compute T/S/U/V for the top boundary layer 
    229127            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 
     
    239137            CALL sbc_isf_cav (kt) 
    240138 
    241          ELSE IF (nn_isf == 2) THEN 
    242             ! Beckmann and Goosse parametrisation  
     139         CASE ( 2 )    ! Beckmann and Goosse parametrisation  
    243140            stbl(:,:)   = soce 
    244141            CALL sbc_isf_bg03(kt) 
    245142 
    246          ELSE IF (nn_isf == 3) THEN 
    247             ! specified runoff in depth (Mathiot et al., XXXX in preparation) 
     143         CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation) 
    248144            CALL fld_read ( kt, nn_fsbc, sf_rnfisf   ) 
    249145            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fwf  flux from the isf (fwfisf <0 mean melting)  
    250             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat flux 
     146            qisf(:,:)   = fwfisf(:,:) * rlfusisf             ! heat flux 
    251147            stbl(:,:)   = soce 
    252148 
    253          ELSE IF (nn_isf == 4) THEN 
    254             ! specified fwf and heat flux forcing beneath the ice shelf 
     149         CASE ( 4 )    ! specified fwf and heat flux forcing beneath the ice shelf 
    255150            CALL fld_read ( kt, nn_fsbc, sf_fwfisf   ) 
    256151            fwfisf(:,:) = - sf_fwfisf(1)%fnow(:,:,1)           ! fwf  flux from the isf (fwfisf <0 mean melting) 
    257             qisf(:,:)   = fwfisf(:,:) * lfusisf                ! heat flux 
     152            qisf(:,:)   = fwfisf(:,:) * rlfusisf               ! heat flux 
    258153            stbl(:,:)   = soce 
    259          END IF 
     154 
     155         END SELECT 
    260156 
    261157         ! compute tsc due to isf 
     
    276172         CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.) 
    277173         CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.) 
    278          CALL lbc_lnk(fwfisf(:,:)   ,'T',1.) 
    279          CALL lbc_lnk(qisf(:,:)     ,'T',1.) 
     174         CALL lbc_lnk(fwfisf(:,:)         ,'T',1.) 
     175         CALL lbc_lnk(qisf(:,:)           ,'T',1.) 
    280176 
    281177         IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    ! 
     
    290186               risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 
    291187            END IF 
    292          ENDIF 
     188         END IF 
    293189         !  
    294190         ! output 
    295          IF( iom_use('qisf'  ) )   CALL iom_put('qisf'  , qisf) 
    296          IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf) 
     191         ! JMM : iom_use not necessary here, qisf and fwfisf are always computed.  
     192         !   If not required in iodef.xml, iom_put does not do anything 
     193!        IF( iom_use('qisf'  ) )   CALL iom_put('qisf'  , qisf) 
     194!        IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf) 
     195         CALL iom_put('qisf'  , qisf) 
     196         CALL iom_put('fwfisf', fwfisf) 
     197 
     198         ! deallocation 
     199         CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  ) 
    297200      END IF 
    298  
    299       ! deallocation 
    300       CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  ) 
    301201   
    302202  END SUBROUTINE sbc_isf 
     
    315215               &    STAT= sbc_isf_alloc ) 
    316216         ! 
    317          IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc ) 
     217         IF( lk_mpp             )   CALL mpp_sum ( sbc_isf_alloc ) 
    318218         IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.') 
    319219         ! 
    320       ENDIF 
     220      END IF 
    321221  END FUNCTION 
    322222 
     223  SUBROUTINE sbc_isf_init 
     224      !!--------------------------------------------------------------------- 
     225      !!                  ***  ROUTINE sbc_isf_init  *** 
     226      !! 
     227      !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation 
     228      !! 
     229      !! ** Method  :  4 parameterizations are available according to nn_isf  
     230      !!               nn_isf = 1 : Realistic ice_shelf formulation 
     231      !!                        2 : Beckmann & Goose parameterization 
     232      !!                        3 : Specified runoff in deptht (Mathiot & al. ) 
     233      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf 
     234      !!---------------------------------------------------------------------- 
     235      INTEGER               :: ji, jj, jk           ! loop index 
     236      INTEGER               :: ik                   ! current level index 
     237      INTEGER               :: ikt, ikb             ! top and bottom level of the isf boundary layer 
     238      INTEGER               :: inum, ierror 
     239      INTEGER               :: ios                  ! Local integer output status for namelist read 
     240      REAL(wp)              :: zhk 
     241      CHARACTER(len=256)    :: cvarzisf, cvarhisf   ! name for isf file 
     242      CHARACTER(LEN=32 )    :: cvarLeff             ! variable name for efficient Length scale 
     243      !!---------------------------------------------------------------------- 
     244      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, & 
     245                         & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 
     246      !!---------------------------------------------------------------------- 
     247 
     248      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs  
     249      READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) 
     250901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp ) 
     251 
     252      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs 
     253      READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 ) 
     254902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 
     255      IF(lwm) WRITE ( numond, namsbc_isf ) 
     256 
     257      IF ( lwp ) WRITE(numout,*) 
     258      IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf' 
     259      IF ( lwp ) WRITE(numout,*) '~~~~~~~~~' 
     260      IF ( lwp ) WRITE(numout,*) 'sbcisf :'  
     261      IF ( lwp ) WRITE(numout,*) '~~~~~~~~' 
     262      IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf 
     263      IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk 
     264      IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl 
     265      IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk  
     266      IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0   
     267      IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0   
     268      IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2  
     269      ! 
     270      ! Allocate public variable 
     271      IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 
     272      ! 
     273      ! initialisation 
     274      qisf(:,:)        = 0._wp  ; fwfisf  (:,:) = 0._wp 
     275      risf_tsc(:,:,:)  = 0._wp  ; fwfisf_b(:,:) = 0._wp 
     276      ! 
     277      ! define isf tbl tickness, top and bottom indice 
     278      SELECT CASE ( nn_isf ) 
     279      CASE ( 1 )  
     280         rhisf_tbl(:,:) = rn_hisf_tbl 
     281         misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv 
     282 
     283      CASE ( 2 , 3 ) 
     284         ALLOCATE( sf_rnfisf(1), STAT=ierror ) 
     285         ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 
     286         CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
     287 
     288         !  read effective lenght (BG03) 
     289         IF (nn_isf == 2) THEN 
     290            CALL iom_open( sn_Leff_isf%clname, inum ) 
     291            cvarLeff = TRIM(sn_Leff_isf%clvar) 
     292            CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 
     293            CALL iom_close(inum) 
     294            ! 
     295            risfLeff = risfLeff*1000.0_wp           !: convertion in m 
     296         END IF 
     297 
     298         ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth) 
     299         CALL iom_open( sn_depmax_isf%clname, inum ) 
     300         cvarhisf = TRIM(sn_depmax_isf%clvar) 
     301         CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base 
     302         CALL iom_close(inum) 
     303         ! 
     304         CALL iom_open( sn_depmin_isf%clname, inum ) 
     305         cvarzisf = TRIM(sn_depmin_isf%clvar) 
     306         CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base 
     307         CALL iom_close(inum) 
     308         ! 
     309         rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer 
     310 
     311         !! compute first level of the top boundary layer 
     312         DO ji = 1, jpi 
     313            DO jj = 1, jpj 
     314                ik = 2 
     315                DO WHILE ( ik <= mbkt(ji,jj) .AND. fsdepw(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO 
     316                misfkt(ji,jj) = ik-1 
     317            END DO 
     318         END DO 
     319 
     320      CASE ( 4 )  
     321         ! as in nn_isf == 1 
     322         rhisf_tbl(:,:) = rn_hisf_tbl 
     323         misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv 
     324          
     325         ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 
     326         ALLOCATE( sf_fwfisf(1), STAT=ierror ) 
     327         ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 
     328         CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
     329 
     330      END SELECT 
     331          
     332      rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 
     333 
     334      ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 
     335      DO jj = 1,jpj 
     336         DO ji = 1,jpi 
     337            ikt = misfkt(ji,jj) 
     338            ikb = misfkt(ji,jj) 
     339            ! thickness of boundary layer at least the top level thickness 
     340            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt)) 
     341 
     342            ! determine the deepest level influenced by the boundary layer 
     343            DO jk = ikt+1, mbkt(ji,jj) 
     344               IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
     345            END DO 
     346            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 
     347            misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl 
     348            r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 
     349 
     350            zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 
     351            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer 
     352         END DO 
     353      END DO 
     354 
     355  END SUBROUTINE sbc_isf_init 
     356 
    323357  SUBROUTINE sbc_isf_bg03(kt) 
    324    !!========================================================================== 
    325    !!                 *** SUBROUTINE sbcisf_bg03  *** 
    326    !! add net heat and fresh water flux from ice shelf melting 
    327    !! into the adjacent ocean using the parameterisation by 
    328    !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 
    329    !!     interaction for climate models", Ocean Modelling 5(2003) 157-170. 
    330    !!  (hereafter BG) 
    331    !!========================================================================== 
    332    !!---------------------------------------------------------------------- 
    333    !!   sbc_isf_bg03      : routine called from sbcmod 
    334    !!---------------------------------------------------------------------- 
    335    !! 
    336    !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting 
    337    !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling 
    338    !! 
    339    !! History : 
    340    !!      !  06-02  (C. Wang) Original code 
    341    !!---------------------------------------------------------------------- 
    342  
    343     INTEGER, INTENT ( in ) :: kt 
    344  
    345     INTEGER :: ji, jj, jk !temporary integer 
    346  
    347     REAL(wp) :: zt_sum    ! sum of the temperature between 200m and 600m 
    348     REAL(wp) :: zt_ave    ! averaged temperature between 200m and 600m 
    349     REAL(wp) :: zt_frz    ! freezing point temperature at depth z 
    350     REAL(wp) :: zpress    ! pressure to compute the freezing point in depth 
    351      
    352     !!---------------------------------------------------------------------- 
    353     IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 
    354      ! 
    355  
    356     ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run ) 
    357     DO ji = 1, jpi 
    358        DO jj = 1, jpj 
    359           jk = misfkt(ji,jj) 
    360           !! Initialize arrays to 0 (each step) 
    361           zt_sum = 0.e0_wp 
    362           IF ( jk .GT. 1 ) THEN 
    363     ! 1. -----------the average temperature between 200m and 600m --------------------- 
    364              DO jk = misfkt(ji,jj),misfkb(ji,jj) 
    365              ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 
    366              ! after verif with UNESCO, wrong sign in BG eq. 2 
    367              ! Calculate freezing temperature 
    368                 CALL eos_fzp(stbl(ji,jj), zt_frz, zpress)  
    369                 zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp 
    370              ENDDO 
    371              zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 
    372      
    373     ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 
    374           ! For those corresponding to zonal boundary     
    375              qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  & 
    376                          & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,jk)  
     358      !!--------------------------------------------------------------------- 
     359      !!                  ***  ROUTINE sbc_isf_bg03  *** 
     360      !! 
     361      !! ** Purpose : add net heat and fresh water flux from ice shelf melting 
     362      !!          into the adjacent ocean 
     363      !! 
     364      !! ** Method  :   See reference 
     365      !! 
     366      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 
     367      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170. 
     368      !!         (hereafter BG) 
     369      !! History : 
     370      !!         06-02  (C. Wang) Original code 
     371      !!---------------------------------------------------------------------- 
     372      INTEGER, INTENT ( in ) :: kt 
     373      ! 
     374      INTEGER  :: ji, jj, jk ! dummy loop index 
     375      INTEGER  :: ik         ! current level 
     376      REAL(wp) :: zt_sum     ! sum of the temperature between 200m and 600m 
     377      REAL(wp) :: zt_ave     ! averaged temperature between 200m and 600m 
     378      REAL(wp) :: zt_frz     ! freezing point temperature at depth z 
     379      REAL(wp) :: zpress     ! pressure to compute the freezing point in depth 
     380      !!---------------------------------------------------------------------- 
     381 
     382      IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 
     383      ! 
     384      DO ji = 1, jpi 
     385         DO jj = 1, jpj 
     386            ik = misfkt(ji,jj) 
     387            !! Initialize arrays to 0 (each step) 
     388            zt_sum = 0.e0_wp 
     389            IF ( ik > 1 ) THEN 
     390               ! 1. -----------the average temperature between 200m and 600m --------------------- 
     391               DO jk = misfkt(ji,jj),misfkb(ji,jj) 
     392                  ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 
     393                  ! after verif with UNESCO, wrong sign in BG eq. 2 
     394                  ! Calculate freezing temperature 
     395                  CALL eos_fzp(stbl(ji,jj), zt_frz, zpress)  
     396                  zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp 
     397               END DO 
     398               zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 
     399               ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 
     400               ! For those corresponding to zonal boundary     
     401               qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  & 
     402                           & * r1_e12t(ji,jj) * tmask(ji,jj,jk) 
    377403              
    378              fwfisf(ji,jj) = qisf(ji,jj) / lfusisf          !fresh water flux kg/(m2s)                   
    379              fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 
    380              !add to salinity trend 
    381           ELSE 
    382              qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 
    383           END IF 
    384        ENDDO 
    385     ENDDO 
    386     ! 
    387     IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03') 
     404               fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !fresh water flux kg/(m2s)                   
     405               fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 
     406               !add to salinity trend 
     407            ELSE 
     408               qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 
     409            END IF 
     410         END DO 
     411      END DO 
     412      ! 
     413      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03') 
     414 
    388415  END SUBROUTINE sbc_isf_bg03 
    389416 
    390    SUBROUTINE sbc_isf_cav( kt ) 
     417  SUBROUTINE sbc_isf_cav( kt ) 
    391418      !!--------------------------------------------------------------------- 
    392419      !!                     ***  ROUTINE sbc_isf_cav  *** 
     
    403430      INTEGER, INTENT(in)          ::   kt         ! ocean time step 
    404431      ! 
    405       REAL(wp), DIMENSION(:,:), POINTER ::   zfrz 
    406       REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas  
    407       REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b 
     432      INTEGER  ::   ji, jj     ! dummy loop indices 
     433      INTEGER  ::   nit 
    408434      REAL(wp) ::   zlamb1, zlamb2, zlamb3 
    409435      REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 
     
    411437      REAL(wp) ::   zeps = 1.e-20_wp         
    412438      REAL(wp) ::   zerr 
    413       INTEGER  ::   ji, jj     ! dummy loop indices 
    414       INTEGER  ::   nit 
     439      REAL(wp), DIMENSION(:,:), POINTER ::   zfrz 
     440      REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas  
     441      REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b 
    415442      LOGICAL  ::   lit 
    416443      !!--------------------------------------------------------------------- 
    417       ! 
    418444      ! coeficient for linearisation of potential tfreez 
    419445      ! Crude approximation for pressure (but commonly used) 
    420       zlamb1=-0.0573_wp 
    421       zlamb2= 0.0832_wp 
    422       zlamb3=-7.53e-08 * grav * rau0 
     446      zlamb1 =-0.0573_wp 
     447      zlamb2 = 0.0832_wp 
     448      zlamb3 =-7.53e-08_wp * grav * rau0 
    423449      IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav') 
    424450      ! 
     
    427453 
    428454      ! initialisation 
    429       zgammat(:,:)=rn_gammat0 ; zgammas (:,:)=rn_gammas0; 
    430       zhtflx (:,:)=0.0_wp     ; zhtflx_b(:,:)=0.0_wp    ; 
    431       zfwflx (:,:)=0.0_wp 
     455      zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0 
     456      zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp     
     457      zfwflx (:,:) = 0.0_wp 
    432458 
    433459      ! compute ice shelf melting 
    434460      nit = 1 ; lit = .TRUE. 
    435461      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine 
    436          IF (nn_isfblk == 1) THEN  
    437             ! ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006) 
     462         SELECT CASE ( nn_isfblk ) 
     463         CASE ( 1 )   ! ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006) 
    438464            ! Calculate freezing temperature 
    439465            CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) ) 
     
    446472               DO ji = 1, jpi 
    447473                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj)) 
    448                   zfwflx(ji,jj) = - zhtflx(ji,jj)/lfusisf 
     474                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf 
    449475               END DO 
    450476            END DO 
     
    454480            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
    455481 
    456          ELSE IF (nn_isfblk == 2 ) THEN 
    457             ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)  
     482         CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 
    458483            ! compute gammat every where (2d) 
    459484            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 
     
    464489               DO ji = 1, jpi 
    465490                  ! compute coeficient to solve the 2nd order equation 
    466                   zeps1=rcp*rau0*zgammat(ji,jj) 
    467                   zeps2=lfusisf*rau0*zgammas(ji,jj) 
    468                   zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps) 
    469                   zeps4=zlamb2+zlamb3*risfdep(ji,jj) 
    470                   zeps6=zeps4-ttbl(ji,jj) 
    471                   zeps7=zeps4-tsurf 
    472                   zaqe=zlamb1 * (zeps1 + zeps3) 
    473                   zaqer=0.5/MIN(zaqe,-zeps) 
    474                   zbqe=zeps1*zeps6+zeps3*zeps7-zeps2 
    475                   zcqe=zeps2*stbl(ji,jj) 
    476                   zdis=zbqe*zbqe-4.0*zaqe*zcqe                
     491                  zeps1 = rcp*rau0*zgammat(ji,jj) 
     492                  zeps2 = rlfusisf*rau0*zgammas(ji,jj) 
     493                  zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps) 
     494                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 
     495                  zeps6 = zeps4-ttbl(ji,jj) 
     496                  zeps7 = zeps4-tsurf 
     497                  zaqe  = zlamb1 * (zeps1 + zeps3) 
     498                  zaqer = 0.5_wp/MIN(zaqe,-zeps) 
     499                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2 
     500                  zcqe  = zeps2*stbl(ji,jj) 
     501                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe                
    477502 
    478503                  ! Presumably zdis can never be negative because gammas is very small compared to gammat 
    479504                  ! compute s freeze 
    480505                  zsfrz=(-zbqe-SQRT(zdis))*zaqer 
    481                   IF ( zsfrz .LT. 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 
     506                  IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 
    482507 
    483508                  ! compute t freeze (eq. 22) 
     
    496521            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
    497522 
    498          ENDIF 
     523         END SELECT 
    499524 
    500525         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration) 
    501          IF ( nn_gammablk .LT. 2 ) THEN ; lit = .FALSE. 
     526         IF ( nn_gammablk < 2 ) THEN ; lit = .FALSE. 
    502527         ELSE                            
    503528            ! check total number of iteration 
    504             IF (nit .GE. 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 
    505             ELSE                   ; nit = nit + 1 
    506             ENDIF 
     529            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 
     530            ELSE                 ; nit = nit + 1 
     531            END IF 
    507532 
    508533            ! compute error between 2 iterations 
    509534            ! if needed save gammat and compute zhtflx_b for next iteration 
    510535            zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 
    511             IF ( zerr .LE. 0.01 ) THEN ; lit = .FALSE. 
    512             ELSE                       ; zhtflx_b(:,:) = zhtflx(:,:) 
    513             ENDIF 
     536            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 
     537            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:) 
     538            END IF 
    514539         END IF 
    515540      END DO 
    516541      ! 
    517       ! output 
    518       IF( iom_use('isfgammat') ) CALL iom_put('isfgammat', zgammat) 
    519       IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas) 
     542      ! output  see JMM comment above 
     543!     IF( iom_use('isfgammat') ) CALL iom_put('isfgammat', zgammat) 
     544!     IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas) 
     545      CALL iom_put('isfgammat', zgammat) 
     546      CALL iom_put('isfgammas', zgammas) 
    520547      !  
    521548      CALL wrk_dealloc( jpi,jpj, zfrz  , zgammat, zgammas  ) 
     
    526553   END SUBROUTINE sbc_isf_cav 
    527554 
    528    SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf ) 
     555   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf ) 
    529556      !!---------------------------------------------------------------------- 
    530557      !! ** Purpose    : compute the coefficient echange for heat flux 
     
    535562      !!                Jenkins et al., 2010, JPO, p2298-2312 
    536563      !!--------------------------------------------------------------------- 
    537       REAL(wp), DIMENSION(:,:), INTENT(out) :: gt, gs 
    538       REAL(wp), DIMENSION(:,:), INTENT(in ) :: zqhisf, zqwisf 
    539  
     564      REAL(wp), DIMENSION(:,:), INTENT(out) :: pgt, pgs 
     565      REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf 
     566      ! 
    540567      INTEGER  :: ikt                         
    541       INTEGER  :: ji,jj                      ! loop index 
     568      INTEGER  :: ji, jj                     ! loop index 
    542569      REAL(wp), DIMENSION(:,:), POINTER :: zustar           ! U, V at T point and friction velocity 
    543570      REAL(wp) :: zdku, zdkv                 ! U, V shear  
     
    551578      REAL(wp) :: zdep 
    552579      REAL(wp) :: zeps = 1.0e-20_wp     
    553       REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant 
    554       REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1) 
     580      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant 
     581      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1) 
    555582      REAL(wp), DIMENSION(2) :: zts, zab 
    556583      !!--------------------------------------------------------------------- 
    557584      CALL wrk_alloc( jpi,jpj, zustar ) 
    558585      ! 
    559       IF      ( nn_gammablk == 0 ) THEN 
    560       !! gamma is constant (specified in namelist) 
    561       !! ISOMIP formulation (Hunter et al, 2006) 
    562          gt(:,:) = rn_gammat0 
    563          gs(:,:) = rn_gammas0 
    564  
    565       ELSE IF ( nn_gammablk == 1 ) THEN 
    566       !! gamma is assume to be proportional to u*  
    567       !! Jenkins et al., 2010, JPO, p2298-2312 
    568       !! Adopted by Asay-Davis et al. (2015) 
    569  
    570       !! compute ustar (eq. 24) 
     586      SELECT CASE ( nn_gammablk ) 
     587      CASE ( 0 ) ! gamma is constant (specified in namelist) 
     588         !! ISOMIP formulation (Hunter et al, 2006) 
     589         pgt(:,:) = rn_gammat0 
     590         pgs(:,:) = rn_gammas0 
     591 
     592      CASE ( 1 ) ! gamma is assume to be proportional to u* 
     593         !! Jenkins et al., 2010, JPO, p2298-2312 
     594         !! Adopted by Asay-Davis et al. (2015) 
     595 
     596         !! compute ustar (eq. 24) 
    571597         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 
    572598 
    573       !! Compute gammats 
    574          gt(:,:) = zustar(:,:) * rn_gammat0 
    575          gs(:,:) = zustar(:,:) * rn_gammas0 
     599         !! Compute gammats 
     600         pgt(:,:) = zustar(:,:) * rn_gammat0 
     601         pgs(:,:) = zustar(:,:) * rn_gammas0 
    576602       
    577       ELSE IF ( nn_gammablk == 2 ) THEN 
    578       !! gamma depends of stability of boundary layer 
    579       !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 
    580       !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 
    581       !! compute ustar 
     603      CASE ( 2 ) ! gamma depends of stability of boundary layer 
     604         !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 
     605         !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 
     606         !! compute ustar 
    582607         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 
    583608 
    584       !! compute Pr and Sc number (can be improved) 
    585          zPr =   13.8 
    586          zSc = 2432.0 
    587  
    588       !! compute gamma mole 
    589          zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0 
    590          zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0 
    591  
    592       !! compute gamma 
     609         !! compute Pr and Sc number (can be improved) 
     610         zPr =   13.8_wp 
     611         zSc = 2432.0_wp 
     612 
     613         !! compute gamma mole 
     614         zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp 
     615         zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp 
     616 
     617         !! compute gamma 
    593618         DO ji=2,jpi 
    594619            DO jj=2,jpj 
     
    596621 
    597622               IF (zustar(ji,jj) == 0._wp) THEN           ! only for kt = 1 I think 
    598                   gt = rn_gammat0 
    599                   gs = rn_gammas0 
     623                  pgt = rn_gammat0 
     624                  pgs = rn_gammas0 
    600625               ELSE 
    601       !! compute Rc number (as done in zdfric.F90) 
    602                   zcoef = 0.5 / fse3w(ji,jj,ikt) 
     626                  !! compute Rc number (as done in zdfric.F90) 
     627                  zcoef = 0.5_wp / fse3w(ji,jj,ikt) 
    603628                  !                                            ! shear of horizontal velocity 
    604629                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  & 
     
    609634                  zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 
    610635 
    611       !! compute bouyancy  
     636                  !! compute bouyancy  
    612637                  zts(jp_tem) = ttbl(ji,jj) 
    613638                  zts(jp_sal) = stbl(ji,jj) 
     
    616641                  CALL eos_rab( zts, zdep, zab ) 
    617642                  ! 
    618       !! compute length scale  
    619                   zbuofdep = grav * ( zab(jp_tem) * zqhisf(ji,jj) - zab(jp_sal) * zqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    620  
    621       !! compute Monin Obukov Length 
     643                  !! compute length scale  
     644                  zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     645 
     646                  !! compute Monin Obukov Length 
    622647                  ! Maximum boundary layer depth 
    623648                  zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) - 0.001_wp 
     
    626651                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 
    627652 
    628       !! compute eta* (stability parameter) 
     653                  !! compute eta* (stability parameter) 
    629654                  zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp))) 
    630655 
    631       !! compute the sublayer thickness 
     656                  !! compute the sublayer thickness 
    632657                  zhnu = 5 * znu / zustar(ji,jj) 
    633658 
    634       !! compute gamma turb 
     659                  !! compute gamma turb 
    635660                  zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) & 
    636661                  &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 
    637662 
    638       !! compute gammats 
    639                   gt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 
    640                   gs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 
     663                  !! compute gammats 
     664                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 
     665                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 
    641666               END IF 
    642667            END DO 
    643668         END DO 
    644          CALL lbc_lnk(gt(:,:),'T',1.) 
    645          CALL lbc_lnk(gs(:,:),'T',1.) 
    646       END IF 
     669         CALL lbc_lnk(pgt(:,:),'T',1.) 
     670         CALL lbc_lnk(pgs(:,:),'T',1.) 
     671      END SELECT 
    647672      CALL wrk_dealloc( jpi,jpj, zustar ) 
    648673 
    649    END SUBROUTINE 
    650  
    651    SUBROUTINE sbc_isf_tbl( varin, varout, cptin ) 
     674   END SUBROUTINE sbc_isf_gammats 
     675 
     676   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin ) 
    652677      !!---------------------------------------------------------------------- 
    653678      !!                  ***  SUBROUTINE sbc_isf_tbl  *** 
     
    656681      !! 
    657682      !!---------------------------------------------------------------------- 
    658       REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin 
    659       REAL(wp), DIMENSION(:,:)  , INTENT(out):: varout 
    660        
    661       CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out 
    662  
     683      REAL(wp), DIMENSION(:,:,:), INTENT( in  ) :: pvarin 
     684      REAL(wp), DIMENSION(:,:)  , INTENT( out ) :: pvarout 
     685      CHARACTER(len=1),           INTENT( in  ) :: cd_ptin ! point of variable in/out 
     686      ! 
    663687      REAL(wp) :: ze3, zhk 
    664688      REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl 
    665689 
    666       INTEGER :: ji,jj,jk                   ! loop index 
    667       INTEGER :: ikt,ikb                    ! top and bottom index of the tbl 
    668   
     690      INTEGER :: ji, jj, jk                  ! loop index 
     691      INTEGER :: ikt, ikb                    ! top and bottom index of the tbl 
     692      !!---------------------------------------------------------------------- 
    669693      ! allocation 
    670694      CALL wrk_alloc( jpi,jpj, zhisf_tbl) 
    671695       
    672696      ! initialisation 
    673       varout(:,:)=0._wp 
     697      pvarout(:,:)=0._wp 
    674698    
    675       ! compute U in the top boundary layer at T- point  
    676       IF (cptin == 'U') THEN 
     699      SELECT CASE ( cd_ptin ) 
     700      CASE ( 'U' ) ! compute U in the top boundary layer at T- point  
    677701         DO jj = 1,jpj 
    678702            DO ji = 1,jpi 
    679703               ikt = miku(ji,jj) ; ikb = miku(ji,jj) 
    680             ! thickness of boundary layer at least the top level thickness 
     704               ! thickness of boundary layer at least the top level thickness 
    681705               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt)) 
    682706 
    683             ! determine the deepest level influenced by the boundary layer 
     707               ! determine the deepest level influenced by the boundary layer 
    684708               DO jk = ikt+1, mbku(ji,jj) 
    685                   IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 
     709                  IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 
    686710               END DO 
    687711               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
    688712 
    689             ! level fully include in the ice shelf boundary layer 
     713               ! level fully include in the ice shelf boundary layer 
    690714               DO jk = ikt, ikb - 1 
    691715                  ze3 = fse3u_n(ji,jj,jk) 
    692                   varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
    693                END DO 
    694  
    695             ! level partially include in ice shelf boundary layer  
     716                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
     717               END DO 
     718 
     719               ! level partially include in ice shelf boundary layer  
    696720               zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
    697                varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 
     721               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
    698722            END DO 
    699723         END DO 
    700724         DO jj = 2,jpj 
    701725            DO ji = 2,jpi 
    702                varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji-1,jj)) 
    703             END DO 
    704          END DO 
    705          CALL lbc_lnk(varout,'T',-1.) 
    706       END IF 
    707  
    708       ! compute V in the top boundary layer at T- point  
    709       IF (cptin == 'V') THEN 
     726               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj)) 
     727            END DO 
     728         END DO 
     729         CALL lbc_lnk(pvarout,'T',-1.) 
     730       
     731      CASE ( 'V' ) ! compute V in the top boundary layer at T- point  
    710732         DO jj = 1,jpj 
    711733            DO ji = 1,jpi 
    712734               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 
    713            ! thickness of boundary layer at least the top level thickness 
     735               ! thickness of boundary layer at least the top level thickness 
    714736               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt)) 
    715737 
    716             ! determine the deepest level influenced by the boundary layer 
     738               ! determine the deepest level influenced by the boundary layer 
    717739               DO jk = ikt+1, mbkv(ji,jj) 
    718                   IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) .LT. zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 
     740                  IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 
    719741               END DO 
    720742               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
    721743 
    722             ! level fully include in the ice shelf boundary layer 
     744               ! level fully include in the ice shelf boundary layer 
    723745               DO jk = ikt, ikb - 1 
    724746                  ze3 = fse3v_n(ji,jj,jk) 
    725                   varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
    726                END DO 
    727  
    728             ! level partially include in ice shelf boundary layer  
     747                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
     748               END DO 
     749 
     750               ! level partially include in ice shelf boundary layer  
    729751               zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
    730                varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 
     752               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
    731753            END DO 
    732754         END DO 
    733755         DO jj = 2,jpj 
    734756            DO ji = 2,jpi 
    735                varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji,jj-1)) 
    736             END DO 
    737          END DO 
    738          CALL lbc_lnk(varout,'T',-1.) 
    739       END IF 
    740  
    741       ! compute T in the top boundary layer at T- point  
    742       IF (cptin == 'T') THEN 
     757               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1)) 
     758            END DO 
     759         END DO 
     760         CALL lbc_lnk(pvarout,'T',-1.) 
     761 
     762      CASE ( 'T' ) ! compute T in the top boundary layer at T- point  
    743763         DO jj = 1,jpj 
    744764            DO ji = 1,jpi 
     
    746766               ikb = misfkb(ji,jj) 
    747767 
    748             ! level fully include in the ice shelf boundary layer 
     768               ! level fully include in the ice shelf boundary layer 
    749769               DO jk = ikt, ikb - 1 
    750770                  ze3 = fse3t_n(ji,jj,jk) 
    751                   varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 
    752                END DO 
    753  
    754             ! level partially include in ice shelf boundary layer  
     771                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 
     772               END DO 
     773 
     774               ! level partially include in ice shelf boundary layer  
    755775               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 
    756                varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 
    757             END DO 
    758          END DO 
    759       END IF 
     776               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
     777            END DO 
     778         END DO 
     779      END SELECT 
    760780 
    761781      ! mask mean tbl value 
    762       varout(:,:) = varout(:,:) * ssmask(:,:) 
     782      pvarout(:,:) = pvarout(:,:) * ssmask(:,:) 
    763783 
    764784      ! deallocation 
     
    780800      !! ** Action  :   phdivn   decreased by the runoff inflow 
    781801      !!---------------------------------------------------------------------- 
    782       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence 
    783       !! 
     802      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence 
     803      !  
    784804      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    785805      INTEGER  ::   ikt, ikb  
     
    790810      zfact   = 0.5_wp 
    791811      ! 
    792       IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water 
     812      IF ( lk_vvl ) THEN     ! need to re compute level distribution of isf fresh water 
    793813         DO jj = 1,jpj 
    794814            DO ji = 1,jpi 
     
    800820               ! determine the deepest level influenced by the boundary layer 
    801821               DO jk = ikt+1, mbkt(ji,jj) 
    802                   IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
     822                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
    803823               END DO 
    804824               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     
    820840               DO jk = ikt, ikb - 1 
    821841                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) & 
    822                     &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 
     842                    &              * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 
    823843               END DO 
    824844               ! level partially include in ice shelf boundary layer  
    825845               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) & 
    826                   &             + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)  
     846                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)  
    827847         END DO 
    828848      END DO 
  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r5621 r5944  
    141141         IF( ln_rnf_tem ) THEN                                       ! use runoffs temperature data 
    142142            rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0 
     143            CALL eos_fzp( sss_m(:,:), ztfrz(:,:) ) 
    143144            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp )             ! if missing data value use SST as runoffs temperature 
    144145               rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 
    145146            END WHERE 
    146147            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg 
    147                ztfrz(:,:) = -1.9 !tfreez( sss_m(:,:) ) !PM to be discuss (trouble if sensitivity study) 
    148                rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * lfusisf * r1_rau0_rcp 
     148               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rlfusisf * r1_rau0_rcp 
    149149            END WHERE 
    150150         ELSE                                                        ! use SST as runoffs temperature 
Note: See TracChangeset for help on using the changeset viewer.