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 8637 for branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC – NEMO

Ignore:
Timestamp:
2017-10-18T19:14:32+02:00 (7 years ago)
Author:
gm
Message:

#1911 (ENHANCE-09): PART I.3 - phasing with updated branch dev_r8183_ICEMODEL revision 8626

Location:
branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC
Files:
1 added
1 deleted
10 edited

Legend:

Unmodified
Added
Removed

  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90

    r8586 r8637  
    1616   !!---------------------------------------------------------------------- 
    1717   USE oce             ! ocean dynamics and tracers variables 
    18    USE ice             ! LIM_3 ice variables 
    19    USE icevar 
    20    USE icectl 
     18   USE ice             ! sea-ice: variables 
     19   USE icevar          ! sea-ice: operations 
     20   USE iceitd          ! sea-ice: rebining 
     21   USE icectl          ! sea-ice: control prints 
    2122   USE phycst          ! physical constant 
    2223   USE eosbn2          ! equation of state 
    … …  
    3940 
    4041   !!---------------------------------------------------------------------- 
    41    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     42   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
    4243   !! $Id: bdyice.F90 8306 2017-07-10 10:18:03Z clem $ 
    4344   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    249250      END DO !jl 
    250251      ! 
     252      ! --- In case categories are out of bounds, do a remapping --- ! 
     253      !     i.e. inputs have not the same ice thickness distribution  
     254      !          (set by rn_himean) than the regional simulation 
     255      IF( jpl > 1 )   CALL ice_itd_reb( kt ) 
    251256      !       
    252257      IF( nn_timing == 1 )   CALL timing_stop('bdy_ice_frs') 
    … …  
    273278      !!------------------------------------------------------------------------------ 
    274279      ! 
    275       IF( nn_timing == 1 ) CALL timing_start('bdy_ice_dyn') 
     280      IF( ln_timing )  CALL timing_start('bdy_ice_dyn') 
    276281      ! 
    277282      DO ib_bdy=1, nb_bdy 
    … …  
    350355      END DO 
    351356      ! 
    352       IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_dyn') 
     357      IF( ln_timing )  CALL timing_stop('bdy_ice_dyn') 
    353358      ! 
    354359    END SUBROUTINE bdy_ice_dyn 

  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

    r8586 r8637  
    4040   USE lib_fortran    ! to use key_nosignedzero 
    4141#if defined key_lim3 
    42    USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b 
     42   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su 
    4343   USE icethd_dh      ! for CALL ice_thd_snwblow 
    4444#endif 
    … …  
    9292   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation 
    9393   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant 
    94    REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice 
     94   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice 
    9595   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant 
    9696   ! 
    … …  
    108108   REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements 
    109109!!gm ref namelist initialize it so remove the setting to false below 
    110    LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm and oce-atm depending on ice concentration (from Lupkes et al. JGR2012) 
     110   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012) 
     111   LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015) 
    111112   ! 
    112    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem) 
     113!!cr   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem) 
     114   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau) 
     115   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens) 
     116   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ce_atm                    ! tansfert coefficient for evaporation   (Q_lat) 
     117   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_zu                      ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme) 
     118   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_zu                      ! air spec. hum.  at wind speed height (needed by Lupkes 2015 bulk scheme) 
     119   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 
    113120 
    114121   INTEGER  ::   nblk           ! choice of the bulk algorithm 
    … …  
    132139      !!             ***  ROUTINE sbc_blk_alloc *** 
    133140      !!------------------------------------------------------------------- 
    134       ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc ) 
     141!!cr      ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc ) 
     142      ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 
     143         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 
    135144      ! 
    136145      IF( lk_mpp             )   CALL mpp_sum ( sbc_blk_alloc ) 
    … …  
    164173         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
    165174         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &  
    166          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12 
     175         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 
    167176      !!--------------------------------------------------------------------- 
    168177      ! 
    … …  
    255264         WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac 
    256265         WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))' 
     266         WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12 
     267         WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15 
    257268         ! 
    258269         WRITE(numout,*) 
    … …  
    361372      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw, zqsb        ! long wave and sensible heat fluxes 
    362373      REAL(wp), DIMENSION(jpi,jpj) ::   zqla, zevap       ! latent heat fluxes and evaporation 
    363       REAL(wp), DIMENSION(jpi,jpj) ::   zCd               ! transfer coefficient for momentum      (tau) 
    364       REAL(wp), DIMENSION(jpi,jpj) ::   zCh               ! transfer coefficient for sensible heat (Q_sens) 
    365       REAL(wp), DIMENSION(jpi,jpj) ::   zCe               ! tansfert coefficient for evaporation   (Q_lat) 
    366374      REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin 
    367       REAL(wp), DIMENSION(jpi,jpj) ::   zt_zu             ! air temperature at wind speed height 
    368       REAL(wp), DIMENSION(jpi,jpj) ::   zq_zu             ! air spec. hum.  at wind speed height 
    369375      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    370376      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
    … …  
    418424      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    419425 
    420  
    421  
    422426      ! ----------------------------------------------------------------------------- ! 
    423427      !     II    Turbulent FLUXES                                                    ! 
    … …  
    435439      ! 
    436440      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2 
    437          &                                              zCd, zCh, zCe, zt_zu, zq_zu, zU_zu ) 
     441         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    438442      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0 
    439          &                                              zCd, zCh, zCe, zt_zu, zq_zu, zU_zu ) 
     443         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    440444      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5 
    441          &                                              zCd, zCh, zCe, zt_zu, zq_zu, zU_zu ) 
     445         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    442446      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF 
    443          &                                              zCd, zCh, zCe, zt_zu, zq_zu, zU_zu ) 
     447         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    444448      CASE DEFAULT 
    445449         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
    … …  
    448452      !                          ! Compute true air density : 
    449453      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 
    450          zrhoa(:,:) = rho_air( zt_zu(:,:)             , zq_zu(:,:)             , sf(jp_slp)%fnow(:,:,1) ) 
     454         zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) ) 
    451455      ELSE                                      ! At zt: 
    452456         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    453457      END IF 
    454458 
    455       Cd_oce(:,:) = zCd(:,:)     ! record value of pure ocean-atm. drag (clem) 
     459!!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef. 
     460!!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef. 
    456461 
    457462      DO jj = 1, jpj             ! tau module, i and j component 
    458463         DO ji = 1, jpi 
    459             zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * zCd(ji,jj)   ! using bulk wind speed 
     464            zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed 
    460465            taum  (ji,jj) = zztmp * wndm  (ji,jj) 
    461466            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
    … …  
    492497      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 
    493498         !! q_air and t_air are given at 10m (wind reference height) 
    494          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*zCe(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) )  ! Evaporation, using bulk wind speed 
    495          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*zCh(:,:)*(zst(:,:) - ztpot(:,:)               )   ! Sensible Heat, using bulk wind speed 
     499         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 
     500         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed 
    496501      ELSE 
    497502         !! q_air and t_air are not given at 10m (wind reference height) 
    498503         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 
    499          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*zCe(:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation ! using bulk wind speed 
    500          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*zCh(:,:)*(zst(:,:) - zt_zu(:,:) )   ! Sensible Heat ! using bulk wind speed 
     504         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 
     505         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed 
    501506      ENDIF 
    502507 
    … …  
    505510 
    506511      IF(ln_ctl) THEN 
    507          CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=zCe , clinfo2=' Ce  : ' ) 
    508          CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zCh , clinfo2=' Ch  : ' ) 
     512         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' ) 
     513         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' ) 
    509514         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 
    510515         CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' ) 
    … …  
    574579      !!--------------------------------------------------------------------- 
    575580      INTEGER  ::   ji, jj    ! dummy loop indices 
    576       REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f      ! relative wind module and components at F-point 
    577       REAL(wp) ::   zwndi_t , zwndj_t                ! relative wind components at T-point 
    578       REAL(wp), DIMENSION(jpi,jpj) ::   zCd, zrhoa   ! transfer coefficient for momentum      (tau) 
     581      REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f   ! relative wind module and components at F-point 
     582      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point 
     583      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa     ! transfer coefficient for momentum      (tau) 
    579584      !!--------------------------------------------------------------------- 
    580585      ! 
    581586      IF( ln_timing )   CALL timing_start('blk_ice_tau') 
    582587      ! 
    583       IF( ln_Cd_L12 ) THEN 
    584          CALL Cdn10_Lupkes2012( zCd )  ! air-ice drag = F(ice concentration) (see Lupkes et al., 2012) 
    585       ELSE 
    586          zCd(:,:) = Cd_ice             ! constant air-ice drag 
    587       ENDIF 
    588  
    589       ! local scalars ( place there for vector optimisation purposes) 
    590       ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
    591       !! 
    592       zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
    593  
    594       !!gm brutal.... 
    595       utau_ice  (:,:) = 0._wp 
    596       vtau_ice  (:,:) = 0._wp 
    597       wndm_ice  (:,:) = 0._wp 
    598       !!gm end 
    599  
    600       ! ----------------------------------------------------------------------------- ! 
    601       !    Wind components and module relative to the moving ocean ( U10m - U_ice )   ! 
    602       ! ----------------------------------------------------------------------------- ! 
     588      ! set transfer coefficients to default sea-ice values 
     589      Cd_atm(:,:) = Cd_ice 
     590      Ch_atm(:,:) = Cd_ice 
     591      Ce_atm(:,:) = Cd_ice 
     592 
     593      wndm_ice(:,:) = 0._wp      !!gm brutal.... 
     594 
     595      ! ------------------------------------------------------------ ! 
     596      !    Wind module relative to the moving ice ( U10m - U_ice )   ! 
     597      ! ------------------------------------------------------------ ! 
    603598      SELECT CASE( cp_ice_msh ) 
    604599      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
    … …  
    606601         DO jj = 2, jpjm1 
    607602            DO ji = 2, jpim1   ! B grid : NO vector opt 
    608                ! ... scalar wind at I-point (fld being at T-point) 
    609                zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   & 
    610                   &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj) 
    611                zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    612                   &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
    613                zwnorm_f = zrhoa(ji,jj) * zCd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    614                ! ... ice stress at I-point 
    615                utau_ice(ji,jj) = zwnorm_f * zwndi_f 
    616                vtau_ice(ji,jj) = zwnorm_f * zwndj_f 
    617603               ! ... scalar wind at T-point (fld being at T-point) 
    618604               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
    … …  
    623609            END DO 
    624610         END DO 
    625          CALL lbc_lnk( utau_ice, 'I', -1. ) 
    626          CALL lbc_lnk( vtau_ice, 'I', -1. ) 
    627611         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    628612         ! 
    629613      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
    630          DO jj = 2, jpj 
    631             DO ji = fs_2, jpi   ! vect. opt. 
     614         DO jj = 2, jpjm1 
     615            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    632616               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
    633617               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
    … …  
    635619            END DO 
    636620         END DO 
     621         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
     622         ! 
     623      END SELECT 
     624 
     625      ! Make ice-atm. drag dependent on ice concentration 
     626      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations 
     627         CALL Cdn10_Lupkes2012( Cd_atm ) 
     628         Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical 
     629      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations 
     630         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )  
     631      ENDIF 
     632 
     633!!      CALL iom_put( "Cd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef. 
     634!!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef. 
     635 
     636      ! local scalars ( place there for vector optimisation purposes) 
     637      ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
     638      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
     639 
     640      !!gm brutal.... 
     641      utau_ice  (:,:) = 0._wp 
     642      vtau_ice  (:,:) = 0._wp 
     643      !!gm end 
     644 
     645      ! ------------------------------------------------------------ ! 
     646      !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
     647      ! ------------------------------------------------------------ ! 
     648      SELECT CASE( cp_ice_msh ) 
     649      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
     650         DO jj = 2, jpjm1 
     651            DO ji = 2, jpim1   ! B grid : NO vector opt 
     652               ! ... scalar wind at I-point (fld being at T-point) 
     653               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   & 
     654                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj) 
     655               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
     656                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
     657               ! ... ice stress at I-point 
     658               zwnorm_f = zrhoa(ji,jj) * Cd_atm(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
     659               utau_ice(ji,jj) = zwnorm_f * zwndi_f 
     660               vtau_ice(ji,jj) = zwnorm_f * zwndj_f 
     661            END DO 
     662         END DO 
     663         CALL lbc_lnk( utau_ice, 'I', -1. ) 
     664         CALL lbc_lnk( vtau_ice, 'I', -1. ) 
     665         ! 
     666      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
    637667         DO jj = 2, jpjm1 
    638668            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    639                utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * zCd(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
     669               utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            & 
    640670                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
    641                vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * zCd(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
     671               vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            & 
    642672                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
    643673            END DO 
    … …  
    645675         CALL lbc_lnk( utau_ice, 'U', -1. ) 
    646676         CALL lbc_lnk( vtau_ice, 'V', -1. ) 
    647          CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    648677         ! 
    649678      END SELECT 
    … …  
    684713      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3) 
    685714      REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa 
    686       REAL(wp), DIMENSION(jpi,jpj)     ::   zCd            ! transfer coefficient for momentum      (tau) 
    687715      !!--------------------------------------------------------------------- 
    688716      ! 
    689717      IF( ln_timing )   CALL timing_start('blk_ice_flx') 
    690718      ! 
    691       IF( ln_Cd_L12 ) THEN 
    692          CALL Cdn10_Lupkes2012( zCd )  ! air-ice drag = F(ice concentration) (see Lupkes et al., 2012) 
    693       ELSE 
    694          zCd(:,:) = Cd_ice             ! constant air-ice drag 
    695       ENDIF 
    696  
    697       ! 
    698       ! local scalars ( place there for vector optimisation purposes) 
     719      ! 
     720      ! local scalars 
    699721      zcoef_dqlw   = 4.0 * 0.95 * Stef 
    700722      zcoef_dqla   = -Ls * 11637800. * (-5897.8) 
    … …  
    724746               ! ----------------------------! 
    725747 
    726                ! ... turbulent heat fluxes 
     748               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau 
    727749               ! Sensible Heat 
    728                z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * zCd(ji,jj) * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
     750               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1)) 
    729751               ! Latent Heat 
    730                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * zCd(ji,jj) * wndm_ice(ji,jj)   & 
    731                   &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
     752               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
     753                  &                ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
    732754               ! Latent heat sensitivity for ice (Dqla/Dt) 
    733755               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
    734                   dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * zCd(ji,jj) * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) ) 
     756                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) / zst2 * EXP(-5897.8 / ptsu(ji,jj,jl)) 
    735757               ELSE 
    736758                  dqla_ice(ji,jj,jl) = 0._wp 
    … …  
    738760 
    739761               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    740                z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * zCd(ji,jj) * wndm_ice(ji,jj) 
     762               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) 
    741763 
    742764               ! ----------------------------! 
    … …  
    935957#if defined key_lim3 
    936958 
    937    SUBROUTINE Cdn10_Lupkes2012( pCd ) 
     959   SUBROUTINE Cdn10_Lupkes2012( Cd ) 
    938960      !!---------------------------------------------------------------------- 
    939961      !!                      ***  ROUTINE  Cdn10_Lupkes2012  *** 
    … …  
    965987      !! 
    966988      !!---------------------------------------------------------------------- 
    967       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pCd   ! air-ice drag coefficient 
     989      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
    968990      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp 
    969991      REAL(wp), PARAMETER ::   znu   = 1._wp 
    … …  
    973995      !!---------------------------------------------------------------------- 
    974996      zcoef = znu + 1._wp / ( 10._wp * zbeta ) 
    975       ! 
     997 
    976998      ! generic drag over a cell partly covered by ice 
    977       !!pCd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag 
    978       !!   &       Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag 
    979       !!   &       zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology 
     999      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag 
     1000      !!   &      Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag 
     1001      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology 
    9801002 
    9811003      ! ice-atm drag 
    982       pCd(:,:) = Cd_ice +  &                                                         ! pure ice drag 
    983          &       zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology 
    984       ! 
     1004      Cd(:,:) = Cd_ice +  &                                                         ! pure ice drag 
     1005         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology 
     1006       
    9851007   END SUBROUTINE Cdn10_Lupkes2012 
     1008 
     1009 
     1010   SUBROUTINE Cdn10_Lupkes2015( Cd, Ch ) 
     1011      !!---------------------------------------------------------------------- 
     1012      !!                      ***  ROUTINE  Cdn10_Lupkes2015  *** 
     1013      !! 
     1014      !! ** pUrpose :    1lternative turbulent transfert coefficients formulation 
     1015      !!                 between sea-ice and atmosphere with distinct momentum  
     1016      !!                 and heat coefficients depending on sea-ice concentration  
     1017      !!                 and atmospheric stability (no meltponds effect for now). 
     1018      !!                 
     1019      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015) 
     1020      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 
     1021      !!                 it considers specific skin and form drags (Andreas et al. 2010) 
     1022      !!                 to compute neutral transfert coefficients for both heat and  
     1023      !!                 momemtum fluxes. Atmospheric stability effect on transfert 
     1024      !!                 coefficient is also taken into account following Louis (1979). 
     1025      !! 
     1026      !! ** References : Lupkes et al. JGR 2015 (theory) 
     1027      !!                 Lupkes et al. ECHAM6 documentation 2015 (implementation) 
     1028      !! 
     1029      !!---------------------------------------------------------------------- 
     1030      ! 
     1031      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
     1032      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch 
     1033      REAL(wp), DIMENSION(jpi,jpj)            ::   zst, zqo_sat, zqi_sat 
     1034      ! 
     1035      ! ECHAM6 constants 
     1036      REAL(wp), PARAMETER ::   z0_skin_ice  = 0.69e-3_wp  ! Eq. 43 [m] 
     1037      REAL(wp), PARAMETER ::   z0_form_ice  = 0.57e-3_wp  ! Eq. 42 [m] 
     1038      REAL(wp), PARAMETER ::   z0_ice       = 1.00e-3_wp  ! Eq. 15 [m] 
     1039      REAL(wp), PARAMETER ::   zce10        = 2.80e-3_wp  ! Eq. 41 
     1040      REAL(wp), PARAMETER ::   zbeta        = 1.1_wp      ! Eq. 41 
     1041      REAL(wp), PARAMETER ::   zc           = 5._wp       ! Eq. 13 
     1042      REAL(wp), PARAMETER ::   zc2          = zc * zc 
     1043      REAL(wp), PARAMETER ::   zam          = 2. * zc     ! Eq. 14 
     1044      REAL(wp), PARAMETER ::   zah          = 3. * zc     ! Eq. 30 
     1045      REAL(wp), PARAMETER ::   z1_alpha     = 1._wp / 0.2_wp  ! Eq. 51 
     1046      REAL(wp), PARAMETER ::   z1_alphaf    = z1_alpha    ! Eq. 56 
     1047      REAL(wp), PARAMETER ::   zbetah       = 1.e-3_wp    ! Eq. 26 
     1048      REAL(wp), PARAMETER ::   zgamma       = 1.25_wp     ! Eq. 26 
     1049      REAL(wp), PARAMETER ::   z1_gamma     = 1._wp / zgamma 
     1050      REAL(wp), PARAMETER ::   r1_3         = 1._wp / 3._wp 
     1051      ! 
     1052      INTEGER  ::   ji, jj         ! dummy loop indices 
     1053      REAL(wp) ::   zthetav_os, zthetav_is, zthetav_zu 
     1054      REAL(wp) ::   zrib_o, zrib_i 
     1055      REAL(wp) ::   zCdn_skin_ice, zCdn_form_ice, zCdn_ice 
     1056      REAL(wp) ::   zChn_skin_ice, zChn_form_ice 
     1057      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw 
     1058      REAL(wp) ::   zCdn_form_tmp 
     1059      !!---------------------------------------------------------------------- 
     1060 
     1061      ! Momentum Neutral Transfert Coefficients (should be a constant) 
     1062      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40 
     1063      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7 
     1064      zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details) 
     1065      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32) 
     1066 
     1067      ! Heat Neutral Transfert Coefficients 
     1068      zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) )   ! Eq. 50 + Eq. 52 (cf Lupkes email for details) 
     1069      
     1070      ! Atmospheric and Surface Variables 
     1071      zst(:,:)     = sst_m(:,:) + rt0                                       ! convert SST from Celcius to Kelvin 
     1072      zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)  , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg] 
     1073      zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] 
     1074      ! 
     1075      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility 
     1076         DO ji = fs_2, fs_jpim1 
     1077            ! Virtual potential temperature [K] 
     1078            zthetav_os = zst(ji,jj)   * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean 
     1079            zthetav_is = tm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice 
     1080            zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu 
     1081             
     1082            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 
     1083            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean 
     1084            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice 
     1085             
     1086            ! Momentum and Heat Neutral Transfert Coefficients 
     1087            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40 
     1088            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53  
     1089                        
     1090            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 
     1091            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 
     1092            z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details) 
     1093            IF( zrib_o <= 0._wp ) THEN 
     1094               zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) )  ! Eq. 10 
     1095               zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26 
     1096                  &             )**zgamma )**z1_gamma 
     1097            ELSE 
     1098               zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12 
     1099               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28 
     1100            ENDIF 
     1101             
     1102            IF( zrib_i <= 0._wp ) THEN 
     1103               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9 
     1104               zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25 
     1105            ELSE 
     1106               zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11 
     1107               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27 
     1108            ENDIF 
     1109             
     1110            ! Momentum Transfert Coefficients (Eq. 38) 
     1111            Cd(ji,jj) = zCdn_skin_ice *   zfmi +  & 
     1112               &        zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 
     1113             
     1114            ! Heat Transfert Coefficients (Eq. 49) 
     1115            Ch(ji,jj) = zChn_skin_ice *   zfhi +  & 
     1116               &        zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 
     1117            ! 
     1118         END DO 
     1119      END DO 
     1120      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. ) 
     1121      ! 
     1122   END SUBROUTINE Cdn10_Lupkes2015 
     1123 
    9861124#endif 
    987     
    9881125 
    9891126   !!====================================================================== 

  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90

    r7646 r8637  
    4444   USE lib_fortran    ! to use key_nosignedzero 
    4545 
    46  
    4746   IMPLICIT NONE 
    4847   PRIVATE 
    … …  
    5049   PUBLIC ::   TURB_COARE   ! called by sbcblk.F90 
    5150 
    52    !! COARE own values for given constants: 
    53    REAL(wp), PARAMETER :: & 
    54       &   zi0     = 600.,  &   !: scale height of the atmospheric boundary layer...1 
    55       &  Beta0    = 1.25, &   !: gustiness parameter 
    56       &  rctv0    = 0.608     !: constant to obtain virtual temperature... 
    57  
    58  
     51   !                                           !! COARE own values for given constants: 
     52   REAL(wp), PARAMETER ::   zi0     = 600.      ! scale height of the atmospheric boundary layer...1 
     53   REAL(wp), PARAMETER ::   Beta0   =   1.25    ! gustiness parameter 
     54   REAL(wp), PARAMETER ::   rctv0   =   0.608   ! constant to obtain virtual temperature... 
     55 
     56   !!---------------------------------------------------------------------- 
    5957CONTAINS 
    6058 
    6159   SUBROUTINE turb_coare( zt, zu, sst, t_zt, ssq, q_zt, U_zu, & 
    62       &                   Cd, Ch, Ce, t_zu, q_zu, U_blk ) 
     60      &                   Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
     61      &                   Cdn, Chn, Cen                       ) 
    6362      !!---------------------------------------------------------------------- 
    6463      !!                      ***  ROUTINE  turb_coare  *** 
    … …  
    106105      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg] 
    107106      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s] 
     107      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    108108      ! 
    109109      INTEGER :: j_itt 
    … …  
    246246      Ce   = ztmp0*q_star/dq_zu 
    247247      ! 
     248      ztmp1 = zu + z0 
     249      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 )) 
     250      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t)) 
     251      Cen = Chn 
     252      ! 
    248253      CALL wrk_dealloc( jpi,jpj, u_star, t_star, q_star, zeta_u, dt_zu, dq_zu ) 
    249254      CALL wrk_dealloc( jpi,jpj, znu_a, z0, z0t, ztmp0, ztmp1, ztmp2 ) 

  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare3p5.F90

    r7646 r8637  
    5858CONTAINS 
    5959 
    60    SUBROUTINE turb_coare3p5( zt, zu, sst, t_zt, ssq, q_zt, U_zu, & 
    61       &                   Cd, Ch, Ce, t_zu, q_zu, U_blk ) 
     60   SUBROUTINE turb_coare3p5( zt, zu, sst, t_zt, ssq, q_zt, U_zu,  & 
     61      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,       & 
     62      &                      Cdn, Chn, Cen                        ) 
    6263      !!---------------------------------------------------------------------------------- 
    6364      !!                      ***  ROUTINE  turb_coare3p5  *** 
    … …  
    105106      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu             [kg/kg] 
    106107      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s] 
     108      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    107109      ! 
    108110      INTEGER :: j_itt 
    … …  
    252254      Ce   = ztmp0*q_star/dq_zu 
    253255      ! 
     256      ztmp1 = zu + z0 
     257      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 )) 
     258      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t)) 
     259      Cen = Chn 
     260      ! 
    254261      CALL wrk_dealloc( jpi,jpj, u_star, t_star, q_star, zeta_u, dt_zu, dq_zu ) 
    255262      CALL wrk_dealloc( jpi,jpj, znu_a, z0, z0t, ztmp0, ztmp1, ztmp2 ) 

  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ecmwf.F90

    r7646 r8637  
    6464 
    6565   SUBROUTINE TURB_ECMWF( zt, zu, sst, t_zt, ssq , q_zt , U_zu,   & 
    66       &                   Cd, Ch, Ce , t_zu, q_zu, U_blk ) 
     66      &                   Cd, Ch, Ce , t_zu, q_zu, U_blk,         & 
     67      &                   Cdn, Chn, Cen                           ) 
    6768      !!---------------------------------------------------------------------------------- 
    6869      !!                      ***  ROUTINE  turb_ecmwf  *** 
    … …  
    112113      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg] 
    113114      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s] 
     115      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    114116      ! 
    115117      INTEGER :: j_itt 
    … …  
    266268            dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 
    267269            dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 
     270 
    268271         END IF 
    269272 
    … …  
    271274         ztmp1 = zu + z0 
    272275         ztmp0 = ztmp1*Linv 
    273          func_m = log(ztmp1) - LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0*Linv) 
     276         func_m = log(ztmp1) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 
    274277         func_h = log(ztmp1) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 
    275278 
    … …  
    280283      ztmp1 = log((zu + z0)/z0q) - psi_h_ecmwf((zu + z0)*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q 
    281284      Ce = vkarmn*vkarmn/(func_m*ztmp1) 
     285 
     286      ztmp1 = zu + z0 
     287      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 )) 
     288      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t)) 
     289      Cen = vkarmn*vkarmn / (log(ztmp1/z0q)*log(ztmp1/z0q)) 
    282290 
    283291      CALL wrk_dealloc( jpi,jpj,   u_star, t_star, q_star, func_m, func_h, dt_zu, dq_zu, Linv ) 

  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ncar.F90

    r7753 r8637  
    5353 
    5454   SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, & 
    55       &                  Cd, Ch, Ce, t_zu, q_zu, U_blk ) 
     55      &                  Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
     56      &                  Cdn, Chn, Cen                       ) 
    5657      !!---------------------------------------------------------------------------------- 
    5758      !!                      ***  ROUTINE  turb_ncar  *** 
    … …  
    112113      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg] 
    113114      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s] 
     115      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    114116      ! 
    115117      INTEGER ::   j_itt 
    … …  
    199201            ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 
    200202            ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10 
     203            Cdn(:,:) = ztmp0 
    201204            sqrt_Cd_n10 = sqrt(ztmp0) 
    202205 
    203206            stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability 
    204207            Cx_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10) 
     208            Chn(:,:) = Cx_n10 
    205209 
    206210            !! Update of transfer coefficients: 
    … …  
    216220 
    217221         Cx_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10 
     222         Cen(:,:) = Cx_n10 
    218223         ztmp1 = 1. + Cx_n10*ztmp0 
    219224         Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    220  
     225         ! 
    221226      END DO 
    222  
     227      ! 
    223228      CALL wrk_dealloc( jpi,jpj,   Cx_n10, sqrt_Cd_n10, zeta_u, stab ) 
    224229      CALL wrk_dealloc( jpi,jpj,   zpsi_h_u, ztmp0, ztmp1, ztmp2 ) 
    225  
     230      ! 
    226231      IF( nn_timing == 1 )   CALL timing_stop('turb_ncar') 
    227  
     232      ! 
    228233   END SUBROUTINE turb_ncar 
    229234 

  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r8586 r8637  
    3333   USE geo2ocean      !  
    3434   USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev 
    35    USE albedooce      !  
     35   USE ocealb         !  
    3636   USE eosbn2         !  
    3737   USE sbcrnf, ONLY : l_rnfcpl 
    … …  
    173173                                         !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel) 
    174174   TYPE ::   DYNARR      
    175       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z3    
     175      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3    
    176176   END TYPE DYNARR 
    177177 
    178178   TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                     ! all fields recieved from the atmosphere 
    179179 
    180    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky) 
     180   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   alb_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky) 
    181181 
    182182   REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2]  
    … …  
    202202      ierr(:) = 0 
    203203      ! 
    204       ALLOCATE( albedo_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) ) 
     204      ALLOCATE( alb_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) ) 
    205205       
    206206#if ! defined key_lim3 && ! defined key_cice 
    … …  
    736736      !     2. receiving mixed oce-ice solar radiation  
    737737      IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN 
    738          CALL albedo_oce( zaos, zacs ) 
     738         CALL oce_alb( zaos, zacs ) 
    739739         ! Due to lack of information on nebulosity : mean clear/overcast sky 
    740          albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5 
     740         alb_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5 
    741741      ENDIF 
    742742 
    … …  
    18851885!       ( see OASIS3 user guide, 5th edition, p39 ) 
    18861886         zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   & 
    1887             &            / (  1.- ( albedo_oce_mix(:,:  ) * ziceld(:,:)       & 
    1888             &                     + palbi         (:,:,1) * picefr(:,:) ) ) 
     1887            &            / (  1.- ( alb_oce_mix(:,:  ) * ziceld(:,:)       & 
     1888            &                     + palbi      (:,:,1) * picefr(:,:) ) ) 
    18891889      END SELECT 
    18901890      IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle 
    … …  
    20522052                   ztmp1(:,:) = SUM( alb_ice (:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 ) / SUM( a_i(:,:,1:jpl), dim=3 ) 
    20532053                ELSEWHERE 
    2054                    ztmp1(:,:) = albedo_oce_mix(:,:) 
     2054                   ztmp1(:,:) = alb_oce_mix(:,:) 
    20552055                END WHERE 
    20562056             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%clcat' ) 
    … …  
    20802080 
    20812081      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean 
    2082          ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:) 
     2082         ztmp1(:,:) = alb_oce_mix(:,:) * zfr_l(:,:) 
    20832083         DO jl=1,jpl 
    20842084            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl) 
    2085          ENDDO 
     2085         END DO 
    20862086         CALL cpl_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
    20872087      ENDIF 

  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90

    r8586 r8637  
    236236      ENDIF 
    237237      ! 
    238       IF( .NOT. l_ssm_mean ) THEN   ! default initialisation. needed by lim_istate 
     238      IF( .NOT.l_ssm_mean ) THEN   ! default initialisation. needed by iceistate 
    239239         ! 
    240240         IF(lwp) WRITE(numout,*) '   default initialisation of ss._m arrays' 

  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90

    r8568 r8637  
    3333   PRIVATE 
    3434 
    35    PUBLIC   zdf_iwm         ! called in step module  
    36    PUBLIC   zdf_iwm_init    ! called in nemogcm module  
    37  
    38    !                       !!* Namelist  namzdf_iwm : internal wave-driven mixing * 
    39    INTEGER  ::  nn_zpyc     ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2) 
    40    LOGICAL  ::  ln_mevar    ! variable (=T) or constant (=F) mixing efficiency 
    41    LOGICAL  ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F) 
    42  
    43    REAL(wp) ::  r1_6 = 1._wp / 6._wp 
    44  
    45    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ebot_iwm     ! power available from high-mode wave breaking (W/m2) 
    46    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   epyc_iwm     ! power available from low-mode, pycnocline-intensified wave breaking (W/m2) 
    47    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ecri_iwm     ! power available from low-mode, critical slope wave breaking (W/m2) 
    48    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbot_iwm     ! WKB decay scale for high-mode energy dissipation (m) 
    49    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcri_iwm     ! decay scale for low-mode critical slope dissipation (m) 
    50    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   emix_iwm     ! local energy density available for mixing (W/kg) 
    51    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bflx_iwm     ! buoyancy flux Kz * N^2 (W/kg) 
    52    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   pcmap_iwm    ! vertically integrated buoyancy flux (W/m2) 
    53    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_ratio    ! S/T diffusivity ratio (only for ln_tsdiff=T) 
    54    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_wave     ! Internal wave-induced diffusivity 
     35   PUBLIC   zdf_iwm        ! called in step module  
     36   PUBLIC   zdf_iwm_init   ! called in nemogcm module  
     37 
     38   !                      !!* Namelist  namzdf_iwm : internal wave-driven mixing * 
     39   INTEGER ::  nn_zpyc     ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2) 
     40   LOGICAL ::  ln_mevar    ! variable (=T) or constant (=F) mixing efficiency 
     41   LOGICAL ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F) 
     42 
     43   REAL(wp)::  r1_6 = 1._wp / 6._wp 
     44 
     45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ebot_iwm   ! power available from high-mode wave breaking (W/m2) 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   epyc_iwm   ! power available from low-mode, pycnocline-intensified wave breaking (W/m2) 
     47   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ecri_iwm   ! power available from low-mode, critical slope wave breaking (W/m2) 
     48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbot_iwm   ! WKB decay scale for high-mode energy dissipation (m) 
     49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hcri_iwm   ! decay scale for low-mode critical slope dissipation (m) 
    5550 
    5651   !! * Substitutions 
    … …  
    6762      !!                ***  FUNCTION zdf_iwm_alloc  *** 
    6863      !!---------------------------------------------------------------------- 
    69       ALLOCATE(     ebot_iwm(jpi,jpj),  epyc_iwm(jpi,jpj),  ecri_iwm(jpi,jpj)    ,   & 
    70       &             hbot_iwm(jpi,jpj),  hcri_iwm(jpi,jpj),  emix_iwm(jpi,jpj,jpk),   & 
    71       &         bflx_iwm(jpi,jpj,jpk), pcmap_iwm(jpi,jpj), zav_ratio(jpi,jpj,jpk),   &  
    72       &         zav_wave(jpi,jpj,jpk), STAT=zdf_iwm_alloc     ) 
     64      ALLOCATE( ebot_iwm(jpi,jpj),  epyc_iwm(jpi,jpj),  ecri_iwm(jpi,jpj) ,     & 
     65      &         hbot_iwm(jpi,jpj),  hcri_iwm(jpi,jpj)                     , STAT=zdf_iwm_alloc ) 
    7366      ! 
    7467      IF( lk_mpp             )   CALL mpp_sum ( zdf_iwm_alloc ) 
    … …  
    8578      !! 
    8679      !! ** Method  : - internal wave-driven vertical mixing is given by: 
    87       !!                  Kz_wave = min(  100 cm2/s, f(  Reb = emix_iwm /( Nu * N^2 )  ) 
    88       !!              where emix_iwm is the 3D space distribution of the wave-breaking  
     80      !!                  Kz_wave = min(  100 cm2/s, f(  Reb = zemx_iwm /( Nu * N^2 )  ) 
     81      !!              where zemx_iwm is the 3D space distribution of the wave-breaking  
    8982      !!              energy and Nu the molecular kinematic viscosity. 
    9083      !!              The function f(Reb) is linear (constant mixing efficiency) 
    9184      !!              if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T. 
    9285      !! 
    93       !!              - Compute emix_iwm, the 3D power density that allows to compute 
     86      !!              - Compute zemx_iwm, the 3D power density that allows to compute 
    9487      !!              Reb and therefrom the wave-induced vertical diffusivity. 
    9588      !!              This is divided into three components: 
    9689      !!                 1. Bottom-intensified low-mode dissipation at critical slopes 
    97       !!                     emix_iwm(z) = ( ecri_iwm / rau0 ) * EXP( -(H-z)/hcri_iwm ) 
     90      !!                     zemx_iwm(z) = ( ecri_iwm / rau0 ) * EXP( -(H-z)/hcri_iwm ) 
    9891      !!                                   / ( 1. - EXP( - H/hcri_iwm ) ) * hcri_iwm 
    9992      !!              where hcri_iwm is the characteristic length scale of the bottom  
    10093      !!              intensification, ecri_iwm a map of available power, and H the ocean depth. 
    10194      !!                 2. Pycnocline-intensified low-mode dissipation 
    102       !!                     emix_iwm(z) = ( epyc_iwm / rau0 ) * ( sqrt(rn2(z))^nn_zpyc ) 
     95      !!                     zemx_iwm(z) = ( epyc_iwm / rau0 ) * ( sqrt(rn2(z))^nn_zpyc ) 
    10396      !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) ) 
    10497      !!              where epyc_iwm is a map of available power, and nn_zpyc 
    … …  
    10699      !!              energy dissipation. 
    107100      !!                 3. WKB-height dependent high mode dissipation 
    108       !!                     emix_iwm(z) = ( ebot_iwm / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm) 
     101      !!                     zemx_iwm(z) = ( ebot_iwm / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm) 
    109102      !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w(z) ) 
    110103      !!              where hbot_iwm is the characteristic length scale of the WKB bottom  
    … …  
    121114      !!                     avs  = avt  +    av_wave * diffusivity_ratio(Reb) 
    122115      !! 
    123       !! ** Action  : - Define emix_iwm used to compute internal wave-induced mixing 
    124       !!              - avt, avs, avm, increased by internal wave-driven mixing     
     116      !! ** Action  : - avt, avs, avm, increased by tide internal wave-driven mixing     
    125117      !! 
    126118      !! References :  de Lavergne et al. 2015, JPO; 2016, in prep. 
    … …  
    132124      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    133125      REAL(wp) ::   zztmp        ! scalar workspace 
    134       REAL(wp), DIMENSION(jpi,jpj)     ::  zfact     ! Used for vertical structure 
    135       REAL(wp), DIMENSION(jpi,jpj)     ::  zhdep     ! Ocean depth 
    136       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwkb      ! WKB-stretched height above bottom 
    137       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zweight   ! Weight for high mode vertical distribution 
    138       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  znu_t     ! Molecular kinematic viscosity (T grid) 
    139       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  znu_w     ! Molecular kinematic viscosity (W grid) 
    140       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zReb      ! Turbulence intensity parameter 
     126      REAL(wp), DIMENSION(jpi,jpj)     ::   zfact       ! Used for vertical structure 
     127      REAL(wp), DIMENSION(jpi,jpj)     ::   zhdep       ! Ocean depth 
     128      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwkb        ! WKB-stretched height above bottom 
     129      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zweight     ! Weight for high mode vertical distribution 
     130      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znu_t       ! Molecular kinematic viscosity (T grid) 
     131      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znu_w       ! Molecular kinematic viscosity (W grid) 
     132      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zReb        ! Turbulence intensity parameter 
     133      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zemx_iwm    ! local energy density available for mixing (W/kg) 
     134      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zav_ratio   ! S/T diffusivity ratio (only for ln_tsdiff=T) 
     135      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zav_wave    ! Internal wave-induced diffusivity 
     136      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z3d  ! 3D workspace used for iom_put  
     137      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   z2d  ! 2D     -      -    -     - 
    141138      !!---------------------------------------------------------------------- 
    142139      ! 
    143140      IF( ln_timing )   CALL timing_start('zdf_iwm') 
    144141      ! 
    145       !                      ! ----------------------------- ! 
    146       !                      !  Internal wave-driven mixing  !  (compute zav_wave) 
    147       !                      ! ----------------------------- ! 
     142      !                       !* Set to zero the 1st and last vertical levels of appropriate variables 
     143      zemx_iwm (:,:,1) = 0._wp   ;   zemx_iwm (:,:,jpk) = 0._wp 
     144      zav_ratio(:,:,1) = 0._wp   ;   zav_ratio(:,:,jpk) = 0._wp 
     145      zav_wave (:,:,1) = 0._wp   ;   zav_wave (:,:,jpk) = 0._wp 
     146      ! 
     147      !                       ! ----------------------------- ! 
     148      !                       !  Internal wave-driven mixing  !  (compute zav_wave) 
     149      !                       ! ----------------------------- ! 
    148150      !                              
    149       !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth, 
     151      !                       !* Critical slope mixing: distribute energy over the time-varying ocean depth, 
    150152      !                                                 using an exponential decay from the seafloor. 
    151153      DO jj = 1, jpj                ! part independent of the level 
    … …  
    156158         END DO 
    157159      END DO 
    158  
     160!!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept_n - sshn 
    159161      DO jk = 2, jpkm1              ! complete with the level-dependent part 
    160          emix_iwm(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w_n(:,:,jk  ) - zhdep(:,:) ) / hcri_iwm(:,:) )                      & 
     162         zemx_iwm(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w_n(:,:,jk  ) - zhdep(:,:) ) / hcri_iwm(:,:) )                      & 
    161163            &                             - EXP( ( gde3w_n(:,:,jk-1) - zhdep(:,:) ) / hcri_iwm(:,:) )  ) * wmask(:,:,jk)   & 
    162164            &                          / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
    … …  
    186188         ! 
    187189         DO jk = 2, jpkm1              ! complete with the level-dependent part 
    188             emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
     190            zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    189191         END DO 
    190192         ! 
    … …  
    203205         ! 
    204206         DO jk = 2, jpkm1              ! complete with the level-dependent part 
    205             emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
     207            zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
    206208         END DO 
    207209         ! 
    … …  
    253255      ! 
    254256      DO jk = 2, jpkm1              ! complete with the level-dependent part 
    255          emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   & 
     257         zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   & 
    256258            &                                / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
    257259!!gm  use of e3t_n just above? 
    … …  
    269271      ! Calculate turbulence intensity parameter Reb 
    270272      DO jk = 2, jpkm1 
    271          zReb(:,:,jk) = emix_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) ) 
     273         zReb(:,:,jk) = zemx_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) ) 
    272274      END DO 
    273275      ! 
    … …  
    349351      !                             !* output internal wave-driven mixing coefficient 
    350352      CALL iom_put( "av_wave", zav_wave ) 
    351                                     !* output useful diagnostics: N^2, Kz * N^2 (bflx_iwm),  
    352                                     !  vertical integral of rau0 * Kz * N^2 (pcmap_iwm), energy density (emix_iwm) 
     353                                    !* output useful diagnostics: Kz*N^2 ,  
     354!!gm Kz*N2 should take into account the ratio avs/avt if it is used.... (see diaar5) 
     355                                    !  vertical integral of rau0 * Kz * N^2 , energy density (zemx_iwm) 
    353356      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 
    354          bflx_iwm(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 
    355          pcmap_iwm(:,:) = 0._wp 
     357         ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) ) 
     358         z3d(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 
     359         z2d(:,:) = 0._wp 
    356360         DO jk = 2, jpkm1 
    357             pcmap_iwm(:,:) = pcmap_iwm(:,:) + e3w_n(:,:,jk) * bflx_iwm(:,:,jk) * wmask(:,:,jk) 
    358          END DO 
    359          pcmap_iwm(:,:) = rau0 * pcmap_iwm(:,:) 
    360          CALL iom_put( "bflx_iwm", bflx_iwm ) 
    361          CALL iom_put( "pcmap_iwm", pcmap_iwm ) 
    362       ENDIF 
    363       CALL iom_put( "bn2", rn2 ) 
    364       CALL iom_put( "emix_iwm", emix_iwm ) 
     361            z2d(:,:) = z2d(:,:) + e3w_n(:,:,jk) * z3d(:,:,jk) * wmask(:,:,jk) 
     362         END DO 
     363         z2d(:,:) = rau0 * z2d(:,:) 
     364         CALL iom_put( "bflx_iwm", z3d ) 
     365         CALL iom_put( "pcmap_iwm", z2d ) 
     366         DEALLOCATE( z2d , z3d ) 
     367      ENDIF 
     368      CALL iom_put( "emix_iwm", zemx_iwm ) 
    365369       
    366370      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk) 
    … …  
    466470      ecri_iwm(:,:) = ecri_iwm(:,:) * ssmask(:,:) 
    467471 
    468       ! Set once for all to zero the first and last vertical levels of appropriate variables 
    469       emix_iwm (:,:, 1 ) = 0._wp 
    470       emix_iwm (:,:,jpk) = 0._wp 
    471       zav_ratio(:,:, 1 ) = 0._wp 
    472       zav_ratio(:,:,jpk) = 0._wp 
    473       zav_wave (:,:, 1 ) = 0._wp 
    474       zav_wave (:,:,jpk) = 0._wp 
    475  
    476472      zbot = glob_sum( e1e2t(:,:) * ebot_iwm(:,:) ) 
    477473      zpyc = glob_sum( e1e2t(:,:) * epyc_iwm(:,:) ) 

  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/step.F90

    r8586 r8637  
    6262      !!                     ***  ROUTINE stp  *** 
    6363      !! 
    64       !! ** Purpose : - Time stepping of OPA (momentum and active tracer eqs.) 
    65       !!              - Time stepping of LIM (dynamic and thermodynamic eqs.) 
    66       !!              - Tme stepping  of TRC (passive tracer eqs.) 
     64      !! ** Purpose : - Time stepping of OPA  (momentum and active tracer eqs.) 
     65      !!              - Time stepping of ESIM (dynamic and thermodynamic eqs.) 
     66      !!              - Time stepping of TRC (passive tracer eqs.) 
    6767      !! 
    6868      !! ** Method  : -1- Update forcings and data 
Note: See TracChangeset for help on using the changeset viewer.