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

ISF: change related to reviewers comments

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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 
Note: See TracChangeset for help on using the changeset viewer.