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 13806 for NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC – NEMO

Ignore:
Timestamp:
2020-11-17T16:58:38+01:00 (4 years ago)
Author:
laurent
Message:

Various improvements and cleaning + keep up with "r13801" of the trunk.

Location:
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbc_ice.F90

    r13655 r13806  
    2020# endif 
    2121# if defined key_cice 
    22    USE ice_domain_size, only: ncat  
     22   USE ice_domain_size, only: ncat 
    2323#endif 
    2424   USE lib_mpp          ! MPP library 
     
    3232# if defined  key_si3 
    3333   LOGICAL         , PUBLIC, PARAMETER ::   lk_si3     = .TRUE.   !: SI3 ice model 
    34    LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .FALSE.  !: no CICE  
     34   LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .FALSE.  !: no CICE 
    3535# endif 
    3636# if defined  key_cice 
     
    4747   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   alb_ice        !: ice albedo                                       [-] 
    4848 
    49    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qml_ice        !: heat available for snow / ice surface melting     [W/m2]  
    50    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice        !: heat conduction flux in the layer below surface   [W/m2]  
     49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qml_ice        !: heat available for snow / ice surface melting     [W/m2] 
     50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice        !: heat conduction flux in the layer below surface   [W/m2] 
    5151   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qtr_ice_top    !: solar flux transmitted below the ice surface      [W/m2] 
    5252 
     
    8787   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iu              !: ice fraction at NEMO U point 
    8888   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iv              !: ice fraction at NEMO V point 
    89     
     89 
    9090   ! variables used in the coupled interface 
    9191   INTEGER , PUBLIC, PARAMETER ::   jpl = ncat 
    92    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice  
    93     
     92   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice 
     93 
    9494   ! already defined in ice.F90 for SI3 
    9595   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i 
     
    9898   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tatm_ice       !: air temperature [K] 
    9999#endif 
    100  
    101    !#LB: REAL(wp), PUBLIC, SAVE ::   pp_cldf = 0.81    !: cloud fraction over sea ice, summer CLIO value   [-] !#LB => moved to sbc_phy.F90 !!! 
    102100 
    103101   !! arrays relating to embedding ice in the ocean 
     
    145143         &                     v_ice(jpi,jpj)        , alb_ice(jpi,jpj,1)    , & 
    146144         &                     emp_ice(jpi,jpj)      , qns_ice(jpi,jpj,1)    , dqns_ice(jpi,jpj,1)   , & 
    147          &                     STAT= ierr(3) )       
     145         &                     STAT= ierr(3) ) 
    148146      IF( ln_cpl )   ALLOCATE( h_i(jpi,jpj,jpl) , h_s(jpi,jpj,jpl) , STAT=ierr(4) ) 
    149147#endif 
     
    168166   LOGICAL         , PUBLIC, PARAMETER ::   lk_si3     = .FALSE.  !: no SI3 ice model 
    169167   LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .FALSE.  !: no CICE ice model 
    170    !#LB: REAL(wp)        , PUBLIC, PARAMETER ::   pp_cldf    = 0.81     !: cloud fraction over sea ice, summer CLIO value   [-] !#LB => moved to sbc_phy.F90 !!!   
    171    INTEGER         , PUBLIC, PARAMETER ::   jpl = 1  
     168 
     169   INTEGER         , PUBLIC, PARAMETER ::   jpl = 1 
    172170   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice                        ! jpi, jpj 
    173171   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_ice, alb_ice, qns_ice, dqns_ice  ! (jpi,jpj,jpl) 
     
    178176   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   topmelt, botmelt 
    179177   ! 
    180    !! arrays related to embedding ice in the ocean.  
    181    !! These arrays need to be declared even if no ice model is required.  
     178   !! arrays related to embedding ice in the ocean. 
     179   !! These arrays need to be declared even if no ice model is required. 
    182180   !! In the no ice model or traditional levitating ice cases they contain only zeros 
    183181   !! --------------------- 
  • NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbc_oce.F90

    r13655 r13806  
    159159   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frq_m     !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-] 
    160160 
    161    !#LB: 
    162161   !!---------------------------------------------------------------------- 
    163162   !!                     Surface atmospheric fields 
     
    165164   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_air_zt       !: specific humidity of air at z=zt [kg/kg]ww 
    166165   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: theta_air_zt   !: potential temperature of air at z=zt [K] 
    167    !#LB. 
    168166 
    169167    
  • NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk.F90

    r13719 r13806  
    107107   ! 
    108108   !#LB: 
    109    LOGICAL  ::   ln_Cx_ice_cst      ! use constant ice-air bulk transfer coefficients (value given in namelist's rn_Cd_i, rn_Ce_i & rn_Ch_i) 
    110    REAL(wp) ::   rn_Cd_i, rn_Ce_i, rn_Ch_i 
    111    LOGICAL  ::   ln_Cx_ice_LU12      ! ice-atm drag = F( ice concentration )                        (Lupkes et al. JGR2012) 
    112    LOGICAL  ::   ln_Cx_ice_LG15      ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 
     109   LOGICAL  ::   ln_Cx_ice_cst             ! use constant air-ice bulk transfer coefficients (value given in namelist's rn_Cd_i, rn_Ce_i & rn_Ch_i) 
     110   REAL(wp) ::   rn_Cd_i, rn_Ce_i, rn_Ch_i ! values for  "    " 
     111   LOGICAL  ::   ln_Cx_ice_LU12            ! air-ice bulk transfer coefficients based on Lupkes et al., 2012) 
     112   LOGICAL  ::   ln_Cx_ice_LG15            ! air-ice bulk transfer coefficients based on Lupkes & Gryanik, 2015) 
    113113   !#LB. 
    114114   ! 
     
    137137   LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 
    138138   LOGICAL  ::   ln_humi_rlh    ! humidity read in files ("sn_humi") is relative humidity     [%] if .true. !LB 
    139    LOGICAL  ::   ln_tpot        !!GS: flag to compute or not potential temperature 
     139   LOGICAL  ::   ln_tair_pot    ! temperature read in files ("sn_tair") is already potential temperature (not absolute) 
    140140   ! 
    141141   INTEGER  ::   nhumi          ! choice of the bulk algorithm 
     
    218218      NAMELIST/namsbc_blk/ ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, ln_ANDREAS, &   ! bulk algorithm 
    219219         &                 rn_zqt, rn_zu, nn_iter_algo, ln_skin_cs, ln_skin_wl,       & 
    220          &                 rn_pfac, rn_efac,                                & 
    221          &                 ln_crt_fbk, rn_stau_a, rn_stau_b,                &   ! current feedback 
    222          &                 ln_humi_sph, ln_humi_dpt, ln_humi_rlh, ln_tpot,  & 
    223          &                 ln_Cx_ice_cst, rn_Cd_i, rn_Ce_i, rn_Ch_i,        & 
    224          &                 ln_Cx_ice_LU12, ln_Cx_ice_LG15,                  & 
    225          &                 cn_dir,                                          & 
    226          &                 sn_wndi, sn_wndj, sn_qsr, sn_qlw ,               &   ! input fields 
    227          &                 sn_tair, sn_humi, sn_prec, sn_snow, sn_slp,      & 
     220         &                 rn_pfac, rn_efac,                                          & 
     221         &                 ln_crt_fbk, rn_stau_a, rn_stau_b,                          &   ! current feedback 
     222         &                 ln_humi_sph, ln_humi_dpt, ln_humi_rlh, ln_tair_pot,        & 
     223         &                 ln_Cx_ice_cst, rn_Cd_i, rn_Ce_i, rn_Ch_i,                  & 
     224         &                 ln_Cx_ice_LU12, ln_Cx_ice_LG15,                            & 
     225         &                 cn_dir,                                                    & 
     226         &                 sn_wndi, sn_wndj, sn_qsr, sn_qlw ,                         &   ! input fields 
     227         &                 sn_tair, sn_humi, sn_prec, sn_snow, sn_slp,                & 
    228228         &                 sn_uoatm, sn_voatm, sn_cc, sn_hpgi, sn_hpgj 
    229229 
     
    567567         !      based on adiabatic lapse-rate (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 
    568568         !      (most reanalysis products provide absolute temp., not potential temp.) 
    569          IF( ln_tpot ) THEN 
     569         IF( ln_tair_pot ) THEN 
     570            ! temperature read into file is already potential temperature, do nothing... 
     571            theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1)             
     572         ELSE 
    570573            ! temperature read into file is ABSOLUTE temperature (that's the case for ECMWF products for example...) 
    571             IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => computing air POTENTIAL temperature out of ABSOLUTE temperature!' 
     574            IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => air temperature converted from ABSOLUTE to POTENTIAL!' 
    572575            theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) + gamma_moist( sf(jp_tair )%fnow(:,:,1), q_air_zt(:,:) ) * rn_zqt 
    573          ELSE 
    574             ! temperature read into file is already potential temperature 
    575             theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) 
    576576         ENDIF 
    577577         ! 
     
    683683      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012) 
    684684      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    685       zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 
    686       zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 
    687       ! ... scalar wind at T-point (not masked) 
    688       wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) 
     685         zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 
     686         zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 
     687         ! ... scalar wind at T-point (not masked) 
     688         wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) 
    689689      END_2D 
    690690#else 
    691691      ! ... scalar wind module at T-point (not masked) 
    692692      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    693       wndm(ji,jj) = SQRT(  pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj)  ) 
     693         wndm(ji,jj) = SQRT(  pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj)  ) 
    694694      END_2D 
    695695#endif 
     
    787787      IF( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp 
    788788         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    789          zztmp = zU_zu(ji,jj) 
    790          wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod 
    791          pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 
    792          psen(ji,jj)   = zztmp * zch_oce(ji,jj) 
    793          pevp(ji,jj)   = zztmp * zce_oce(ji,jj) 
    794          rhoa(ji,jj)   = rho_air( ptair(ji,jj), pqair(ji,jj), pslp(ji,jj) ) 
     789            zztmp = zU_zu(ji,jj) 
     790            wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod 
     791            pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 
     792            psen(ji,jj)   = zztmp * zch_oce(ji,jj) 
     793            pevp(ji,jj)   = zztmp * zce_oce(ji,jj) 
     794            rhoa(ji,jj)   = rho_air( ptair(ji,jj), pqair(ji,jj), pslp(ji,jj) ) 
    795795         END_2D 
    796796      ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation 
     
    825825            zstmax = MIN( rn_stau_a * 3._wp + rn_stau_b, 0._wp )   ! set the max value of Stau corresponding to a wind of 3 m/s (<0) 
    826826            DO_2D( 0, 1, 0, 1 )   ! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop 
    827             zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax )   ! stau (<0) must be smaller than zstmax 
    828             ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj  ) + pu(ji,jj) ) - puatm(ji,jj) ) 
    829             ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji  ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) ) 
    830             taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) ) 
     827               zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax )   ! stau (<0) must be smaller than zstmax 
     828               ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj  ) + pu(ji,jj) ) - puatm(ji,jj) ) 
     829               ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji  ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) ) 
     830               taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) ) 
    831831            END_2D 
    832832         ENDIF 
     
    836836         !     Note that coastal wind stress is not used in the code... so this extra care has no effect 
    837837         DO_2D( 0, 0, 0, 0 )              ! start loop at 2, in case ln_crt_fbk = T 
    838          utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj  ) ) & 
    839             &              * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
    840          vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji  ,jj+1) ) & 
    841             &              * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
     838            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj  ) ) & 
     839               &              * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
     840            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji  ,jj+1) ) & 
     841               &              * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
    842842         END_2D 
    843  
    844843 
    845844         IF( ln_crt_fbk ) THEN 
     
    870869   END SUBROUTINE blk_oce_1 
    871870 
    872     
     871 
    873872   SUBROUTINE blk_oce_2( ptair, pdqlw, pprec, psnow, &   ! <<= in 
    874873      &                   ptsk, psen, plat, pevp     )   ! <<= in 
     
    13391338         DO jl = 1, jpl 
    13401339            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    1341             zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness 
    1342             IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 
     1340               zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness 
     1341               IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 
    13431342            END_2D 
    13441343         END DO 
     
    13541353      DO jl = 1, jpl 
    13551354         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    1356          ! 
    1357          zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness 
    1358             &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 
    1359          ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature 
    1360          ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature 
    1361          zqa0    = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux 
    1362          ! 
    1363          DO iter = 1, nit     ! --- Iterative loop 
    1364             zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards) 
    1365             zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget 
    1366             ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update 
    1367          END DO 
    1368          ! 
    1369          ptsu   (ji,jj,jl) = MIN( rt0, ztsu ) 
    1370          qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 
    1371          qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 
    1372          qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) )  & 
    1373             &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 
    1374  
    1375          ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 
    1376          hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) 
     1355            ! 
     1356            zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness 
     1357               &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 
     1358            ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature 
     1359            ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature 
     1360            zqa0    = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux 
     1361            ! 
     1362            DO iter = 1, nit     ! --- Iterative loop 
     1363               zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards) 
     1364               zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget 
     1365               ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update 
     1366            END DO 
     1367            ! 
     1368            ptsu   (ji,jj,jl) = MIN( rt0, ztsu ) 
     1369            qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 
     1370            qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 
     1371            qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) )  & 
     1372               &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 
     1373 
     1374            ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 
     1375            hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) 
    13771376 
    13781377         END_2D 
  • NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_skin_ecmwf.F90

    r13719 r13806  
    160160      REAL(wp) :: zalfa     !: thermal expansion coefficient of sea-water [1/K] 
    161161      REAL(wp) :: zdTwl_b, zdTwl_n  !: temp. diff. between "almost surface (right below viscous layer) and bottom of WL 
    162       REAL(wp) :: zfr, zeta  
    163       REAL(wp) :: zusw, zusw2  
    164       REAL(wp) :: zLa, zfLa  
    165       REAL(wp) :: flg, zwf, zQabs  
     162      REAL(wp) :: zfr, zeta 
     163      REAL(wp) :: zusw, zusw2 
     164      REAL(wp) :: zLa, zfLa 
     165      REAL(wp) :: flg, zwf, zQabs 
    166166      REAL(wp) :: ZA, ZB, zL1, zL2 
    167167      REAL(wp) :: zcst0, zcst1, zcst2, zcst3 
  • NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcmod.F90

    r13655 r13806  
    4646   USE sbcssr         ! surface boundary condition: sea surface restoring 
    4747   USE sbcrnf         ! surface boundary condition: runoffs 
    48    USE sbcapr         ! surface boundary condition: atmo pressure  
     48   USE sbcapr         ! surface boundary condition: atmo pressure 
    4949   USE sbcfwb         ! surface boundary condition: freshwater budget 
    5050   USE icbstp         ! Icebergs 
     
    137137         WRITE(numout,*) '         ocean-atmosphere coupled formulation       ln_cpl        = ', ln_cpl 
    138138         WRITE(numout,*) '         mixed forced-coupled     formulation       ln_mixcpl     = ', ln_mixcpl 
    139 !!gm  lk_oasis is controlled by key_oasis3  ===>>>  It shoud be removed from the namelist  
     139!!gm  lk_oasis is controlled by key_oasis3  ===>>>  It shoud be removed from the namelist 
    140140         WRITE(numout,*) '         OASIS coupling (with atm or sas)           lk_oasis      = ', lk_oasis 
    141141         WRITE(numout,*) '         components of your executable              nn_components = ', nn_components 
     
    162162      IF( .NOT.ln_wave ) THEN 
    163163         ln_sdw = .false. ; ln_cdgw = .false. ; ln_tauwoc = .false. ; ln_tauw = .false. ; ln_stcor = .false. 
    164       ENDIF  
     164      ENDIF 
    165165      IF( ln_sdw ) THEN 
    166166         IF( .NOT.(nn_sdrift==jp_breivik_2014 .OR. nn_sdrift==jp_li_2017 .OR. nn_sdrift==jp_peakfr) ) & 
     
    187187      !                       !**  check option consistency 
    188188      ! 
    189       IF(lwp) WRITE(numout,*)       !* Single / Multi - executable (NEMO / OPA+SAS)  
     189      IF(lwp) WRITE(numout,*)       !* Single / Multi - executable (NEMO / OPA+SAS) 
    190190      SELECT CASE( nn_components ) 
    191191      CASE( jp_iam_nemo ) 
     
    219219      SELECT CASE( nn_ice ) 
    220220      CASE( 0 )                        !- no ice in the domain 
    221       CASE( 1 )                        !- Ice-cover climatology ("Ice-if" model)   
     221      CASE( 1 )                        !- Ice-cover climatology ("Ice-if" model) 
    222222      CASE( 2 )                        !- SI3  ice model 
    223223         IF( .NOT.( ln_blk .OR. ln_cpl .OR. ln_abl .OR. ln_usr ) )   & 
     
    227227            &                   CALL ctl_stop( 'sbc_init : CICE sea-ice model requires ln_blk or ln_cpl or ln_abl or ln_usr = T' ) 
    228228         IF( lk_agrif                                )   & 
    229             &                   CALL ctl_stop( 'sbc_init : CICE sea-ice model not currently available with AGRIF' )  
     229            &                   CALL ctl_stop( 'sbc_init : CICE sea-ice model not currently available with AGRIF' ) 
    230230      CASE DEFAULT                     !- not supported 
    231231      END SELECT 
    232       IF( ln_diurnal .AND. .NOT. ln_blk )   CALL ctl_stop( "sbc_init: diurnal flux processing only implemented for bulk forcing" ) 
     232      IF( ln_diurnal .AND. .NOT. (ln_blk.OR.ln_abl) )   CALL ctl_stop( "sbc_init: diurnal flux processing only implemented for bulk forcing" ) 
    233233      ! 
    234234      !                       !**  allocate and set required variables 
     
    242242      ! 
    243243      IF( sbc_ssr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_ssr arrays' ) 
    244       IF( .NOT.ln_ssr ) THEN               !* Initialize qrp and erp if no restoring  
     244      IF( .NOT.ln_ssr ) THEN               !* Initialize qrp and erp if no restoring 
    245245         qrp(:,:) = 0._wp 
    246246         erp(:,:) = 0._wp 
     
    331331         &   CALL ctl_warn( 'sbc_init : diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' ) 
    332332      ! 
    333     
     333 
    334334      !                       !**  associated modules : initialization 
    335335      ! 
     
    395395      ! 
    396396      REAL(wp) ::     zthscl        ! wd  tanh scale 
    397       REAL(wp), DIMENSION(jpi,jpj) ::  zwdht, zwght  ! wd dep over wd limit, wgt   
     397      REAL(wp), DIMENSION(jpi,jpj) ::  zwdht, zwght  ! wd dep over wd limit, wgt 
    398398 
    399399      !!--------------------------------------------------------------------- 
     
    427427      ! 
    428428      !                                            !==  sbc formulation  ==! 
    429       !                                                    
     429      ! 
    430430      SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition 
    431431      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx) 
    432       CASE( jp_usr   )     ;   CALL usrdef_sbc_oce( kt, Kbb )                        ! user defined formulation  
     432      CASE( jp_usr   )     ;   CALL usrdef_sbc_oce( kt, Kbb )                        ! user defined formulation 
    433433      CASE( jp_flx     )   ;   CALL sbc_flx       ( kt )                             ! flux formulation 
    434434      CASE( jp_blk     ) 
     
    447447      IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! forced-coupled mixed formulation after forcing 
    448448      ! 
    449       IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )              ! Wind stress provided by waves  
     449      IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )              ! Wind stress provided by waves 
    450450      ! 
    451451      !                                            !==  Misc. Options  ==! 
     
    461461      IF( ln_icebergs    )   THEN 
    462462                                     CALL icb_stp( kt )           ! compute icebergs 
    463          ! Icebergs do not melt over the haloes.  
    464          ! So emp values over the haloes are no more consistent with the inner domain values.  
     463         ! Icebergs do not melt over the haloes. 
     464         ! So emp values over the haloes are no more consistent with the inner domain values. 
    465465         ! A lbc_lnk is therefore needed to ensure reproducibility and restartability. 
    466466         ! see ticket #2113 for discussion about this lbc_lnk. 
     
    476476      ! Special treatment of freshwater fluxes over closed seas in the model domain 
    477477      ! Should not be run if ln_diurnal_only 
    478       IF( l_sbc_clo      )   CALL sbc_clo( kt )    
     478      IF( l_sbc_clo      )   CALL sbc_clo( kt ) 
    479479 
    480480!!$!RBbug do not understand why see ticket 667 
     
    482482!!$      CALL lbc_lnk( 'sbcmod', emp, 'T', 1.0_wp ) 
    483483      IF( ll_wd ) THEN     ! If near WAD point limit the flux for now 
    484          zthscl = atanh(rn_wd_sbcfra)                     ! taper frac default is .999  
     484         zthscl = atanh(rn_wd_sbcfra)                     ! taper frac default is .999 
    485485         zwdht(:,:) = ssh(:,:,Kmm) + ht_0(:,:) - rn_wdmin1   ! do this calc of water 
    486486                                                     ! depth above wd limit once 
Note: See TracChangeset for help on using the changeset viewer.