Changeset 12565


Ignore:
Timestamp:
2020-03-17T15:51:51+01:00 (7 months ago)
Author:
smasson
Message:

dev_r12472_ASINTER-05_Masson_CurrentFeedback: code inmplemenation, see #2156

Location:
NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg

    r12495 r12565  
    9494   !                    !  bulk algorithm : 
    9595   ln_NCAR      = .true.    ! "NCAR"      algorithm   (Large and Yeager 2008) 
    96    ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    97    ln_COARE_3p6 = .false.   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
    98    ln_ECMWF     = .false.   ! "ECMWF"     algorithm   (IFS cycle 31) 
    99       ! 
    100       rn_zqt      = 10.       !  Air temperature & humidity reference height (m) 
    101       rn_zu       = 10.       !  Wind vector reference height (m) 
    102       ln_Cd_L12   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2012) 
    103       ln_Cd_L15   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2015) 
    104       rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow) 
    105       rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.) 
    106       rn_vfac     = 0.        !  multiplicative factor for ocean & ice velocity used to 
    107       !                       !  calculate the wind stress (0.=absolute or 1.=relative winds) 
    108       ln_skin_cs = .false.  !  use the cool-skin parameterization (only available in ECMWF and COARE algorithms) !LB 
    109       ln_skin_wl = .false.  !  use the warm-layer        "               "                    " 
    110       ! 
    111       ln_humi_sph = .true.     !  humidity specified below in "sn_humi" is specific humidity     [kg/kg] if .true. 
    112       ln_humi_dpt = .false.    !  humidity specified below in "sn_humi" is dew-point temperature   [K]   if .true. 
    113       ln_humi_rlh = .false.    !  humidity specified below in "sn_humi" is relative humidity       [%]   if .true. 
    11496   ! 
    11597   cn_dir = './'  !  root directory for the bulk data location 
  • NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg

    r12495 r12565  
    9090   !                    !  bulk algorithm : 
    9191   ln_NCAR      = .true.    ! "NCAR"      algorithm   (Large and Yeager 2008) 
    92    ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    93    ln_COARE_3p6 = .false.   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
    94    ln_ECMWF     = .false.   ! "ECMWF"     algorithm   (IFS cycle 31) 
    95       ! 
    96       rn_zqt      = 10.       !  Air temperature & humidity reference height (m) 
    97       rn_zu       = 10.       !  Wind vector reference height (m) 
    98       ln_Cd_L12   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2012) 
    99       ln_Cd_L15   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2015) 
    100       rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow) 
    101       rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.) 
    102       rn_vfac     = 0.        !  multiplicative factor for ocean & ice velocity used to 
    103       !                       !  calculate the wind stress (0.=absolute or 1.=relative winds) 
    104       ln_skin_cs = .false.  !  use the cool-skin parameterization (only available in ECMWF and COARE algorithms) !LB 
    105       ln_skin_wl = .false.  !  use the warm-layer        "               "                    " 
    106       ! 
    107       ln_humi_sph = .true.     !  humidity specified below in "sn_humi" is specific humidity     [kg/kg] if .true. 
    108       ln_humi_dpt = .false.    !  humidity specified below in "sn_humi" is dew-point temperature   [K]   if .true. 
    109       ln_humi_rlh = .false.    !  humidity specified below in "sn_humi" is relative humidity       [%]   if .true. 
    11092   ! 
    11193   cn_dir      = './'      !  root directory for the bulk data location 
  • NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg

    r12495 r12565  
    9494   !                    !  bulk algorithm : 
    9595   ln_NCAR      = .true.    ! "NCAR"      algorithm   (Large and Yeager 2008) 
    96    ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    97    ln_COARE_3p6 = .false.   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
    98    ln_ECMWF     = .false.   ! "ECMWF"     algorithm   (IFS cycle 31) 
    99       ! 
    100       rn_zqt      = 10.       !  Air temperature & humidity reference height (m) 
    101       rn_zu       = 10.       !  Wind vector reference height (m) 
    102       ln_Cd_L12   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2012) 
    103       ln_Cd_L15   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2015) 
    104       rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow) 
    105       rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.) 
    106       rn_vfac     = 0.        !  multiplicative factor for ocean & ice velocity used to 
    107       !                       !  calculate the wind stress (0.=absolute or 1.=relative winds) 
    108       ln_skin_cs = .false.  !  use the cool-skin parameterization (only available in ECMWF and COARE algorithms) !LB 
    109       ln_skin_wl = .false.  !  use the warm-layer        "               "                    " 
    110       ! 
    111       ln_humi_sph = .true.     !  humidity specified below in "sn_humi" is specific humidity     [kg/kg] if .true. 
    112       ln_humi_dpt = .false.    !  humidity specified below in "sn_humi" is dew-point temperature   [K]   if .true. 
    113       ln_humi_rlh = .false.    !  humidity specified below in "sn_humi" is relative humidity       [%]   if .true. 
    11496   ! 
    11597   cn_dir = './'  !  root directory for the bulk data location 
  • NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/cfgs/ORCA2_ICE_ABL/EXPREF/namelist_cfg

    r12495 r12565  
    109109   !                    !  bulk algorithm : 
    110110   ln_NCAR      = .true.    ! "NCAR"      algorithm   (Large and Yeager 2008) 
    111    ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    112    ln_COARE_3p6 = .false.   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
    113    ln_ECMWF     = .false.   ! "ECMWF"     algorithm   (IFS cycle 31) 
    114       rn_zqt      = 10.     !  Air temperature & humidity reference height (m) 
    115       rn_zu       = 10.     !  Wind vector reference height (m) 
    116       ! 
    117       ! Skin is ONLY available in ECMWF and COARE algorithms: 
    118       ln_skin_cs = .false.  !  use the cool-skin parameterization => set nn_fsbc=1 and ln_dm2dc=.true.! 
    119       ln_skin_wl = .false.  !  use the warm-layer        "        => set nn_fsbc=1 and ln_dm2dc=.true.! 
    120       ! 
    121       ln_humi_sph = .true.  !  humidity specified below in "sn_humi" is specific humidity     [kg/kg] if .true. 
    122       ln_humi_dpt = .false. !  humidity specified below in "sn_humi" is dew-point temperature   [K]   if .true. 
    123       ln_humi_rlh = .false. !  humidity specified below in "sn_humi" is relative humidity       [%]   if .true. 
    124111   ! 
    125112   cn_dir = './'  !  root directory for the bulk data location 
     
    131118   sn_tair     = 'tair_drwnlnd_ERAI_L25Z10_GLOBAL_F128R_ana1d',  24., 'tair'    , .false. , .false. , 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bilinear' , ''    , '' 
    132119   sn_humi     = 'humi_drwnlnd_ERAI_L25Z10_GLOBAL_F128R_ana1d',  24., 'humi'    , .false. , .false. , 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bilinear' , ''    , '' 
    133    sn_hpgi     = 'uhpg_drwnlnd_ERAI_L25Z10_GLOBAL_F128R_ana1d',  24., 'uhpg'    , .false. , .false. , 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic'  , 'UG'  , '' 
    134    sn_hpgj     = 'vhpg_drwnlnd_ERAI_L25Z10_GLOBAL_F128R_ana1d',  24., 'vhpg'    , .false. , .false. , 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic'  , 'VG'  , '' 
    135  
    136120   sn_qsr      = 'ncar_rad.15JUNE2009_fill'                    , 24., 'SWDN_MOD', .false. , .true.  ,  'yearly' , 'weights_core_orca2_bilinear_noc.nc'      , ''    , '' 
    137121   sn_qlw      = 'ncar_rad.15JUNE2009_fill'                    , 24., 'LWDN_MOD', .false. , .true.  ,  'yearly' , 'weights_core_orca2_bilinear_noc.nc'      , ''    , '' 
     
    139123   sn_snow     = 'ncar_precip.15JUNE2009_fill'                 , -1., 'SNOW'    , .false. , .true.  ,  'yearly' , 'weights_core_orca2_bilinear_noc.nc'      , ''    , '' 
    140124   sn_slp      = 'slp.15JUNE2009_fill'                         ,  6., 'SLP'     , .false. , .true.  ,  'yearly' , 'weights_core_orca2_bilinear_noc.nc'      , ''    , '' 
     125   sn_hpgi     = 'uhpg_drwnlnd_ERAI_L25Z10_GLOBAL_F128R_ana1d',  24., 'uhpg'    , .false. , .false. , 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic'  , 'UG'  , '' 
     126   sn_hpgj     = 'vhpg_drwnlnd_ERAI_L25Z10_GLOBAL_F128R_ana1d',  24., 'vhpg'    , .false. , .false. , 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic'  , 'VG'  , '' 
    141127/ 
    142128 
  • NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/cfgs/ORCA2_SAS_ICE/EXPREF/namelist_cfg

    r12377 r12565  
    6666   !                    !  bulk algorithm : 
    6767   ln_NCAR      = .true.    ! "NCAR"      algorithm   (Large and Yeager 2008) 
    68    ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    69    ln_COARE_3p6 = .false.   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
    70    ln_ECMWF     = .false.   ! "ECMWF"     algorithm   (IFS cycle 31) 
    71       ! 
    72       rn_zqt      = 10.       !  Air temperature & humidity reference height (m) 
    73       rn_zu       = 10.       !  Wind vector reference height (m) 
    74       ln_Cd_L12   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2012) 
    75       ln_Cd_L15   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2015) 
    76       rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow) 
    77       rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.) 
    78       rn_vfac     = 0.        !  multiplicative factor for ocean & ice velocity used to 
    79       !                       !  calculate the wind stress (0.=absolute or 1.=relative winds) 
    80       ln_skin_cs = .false.  !  use the cool-skin parameterization (only available in ECMWF and COARE algorithms) !LB 
    81       ln_skin_wl = .false.  !  use the warm-layer        "               "                    " 
    82       ! 
    83       ln_humi_sph = .true.     !  humidity specified below in "sn_humi" is specific humidity     [kg/kg] if .true. 
    84       ln_humi_dpt = .false.    !  humidity specified below in "sn_humi" is dew-point temperature   [K]   if .true. 
    85       ln_humi_rlh = .false.    !  humidity specified below in "sn_humi" is relative humidity       [%]   if .true. 
    8668   ! 
    8769   cn_dir      = './'      !  root directory for the bulk data location 
  • NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/cfgs/SHARED/namelist_ref

    r12551 r12565  
    268268      ln_Cd_L12  = .false.  !  air-ice drags = F(ice conc.) (Lupkes et al. 2012) 
    269269      ln_Cd_L15  = .false.  !  air-ice drags = F(ice conc.) (Lupkes et al. 2015) 
    270       !                     !  - module of the mean stress" data 
     270      ln_crt_fbk = .false.     !  Add surface current feedback to the wind stress (Renault et al. 2020, doi: 10.1029/2019MS001715) 
     271         rn_stau_a = -2.9e-3   !     Alpha from eq. 10: Stau = Alpha * Wnd + Beta 
     272         rn_stau_b =  8.0e-3   !     Beta  
    271273      rn_pfac    = 1.       !  multipl. factor for precipitation (total & snow) 
    272274      rn_efac    = 1.       !  multipl. factor for evaporation (0. or 1.) 
    273       rn_vfac    = 0.       !  multipl. factor for ocean & ice velocity 
    274       !                     !  used to calculate the wind stress 
    275       !                     ! (0. => absolute or 1. => relative winds) 
    276275      ln_skin_cs = .false.  !  use the cool-skin parameterization 
    277276      ln_skin_wl = .false.  !  use the warm-layer parameterization 
     
    291290   sn_tair     = 't_10.15JUNE2009_fill'       ,    6.        , 'T_10_MOD',   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    292291   sn_humi     = 'q_10.15JUNE2009_fill'       ,    6.        , 'Q_10_MOD',   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    293    sn_hpgi     = 'NOT USED'                   ,   24.        , 'uhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'UG'     , '' 
    294    sn_hpgj     = 'NOT USED'                   ,   24.        , 'vhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'VG'     , '' 
    295292   sn_prec     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'PRC_MOD1',   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    296293   sn_snow     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'SNOW'    ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    297294   sn_slp      = 'slp.15JUNE2009_fill'        ,    6.        , 'SLP'     ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
     295   sn_uoatm    = 'NOT USED'                   ,    6.        , 'UOATM'   ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , 'Uoceatm', '' 
     296   sn_voatm    = 'NOT USED'                   ,    6.        , 'VOATM'   ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , 'Voceatm', '' 
     297   sn_hpgi     = 'NOT USED'                   ,   24.        , 'uhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'UG'     , '' 
     298   sn_hpgj     = 'NOT USED'                   ,   24.        , 'vhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'VG'     , '' 
    298299/ 
    299300!----------------------------------------------------------------------- 
  • NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/cfgs/WED025/EXPREF/namelist_cfg

    r12495 r12565  
    139139!----------------------------------------------------------------------- 
    140140   !                    !  bulk algorithm : 
    141    ln_NCAR     = .true.   ! "NCAR"      algorithm   (Large and Yeager 2008) 
    142    ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    143    ln_COARE_3p5 = .false.   ! "COARE 3.5" algorithm   (Edson et al. 2013) 
    144    ln_ECMWF    = .false.   ! "ECMWF"     algorithm   (IFS cycle 31) 
    145  
     141   ln_NCAR     = .true. 
     142   ! 
    146143   cn_dir      = './'      !  root directory for the bulk data location 
    147144   !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________! 
  • NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/src/ABL/ablmod.F90

    r12495 r12565  
    1818   USE dom_oce, ONLY  : tmask   
    1919   USE sbc_oce, ONLY  : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1, rhoa 
    20    USE sbcblk         ! use rn_?fac 
     20   USE sbcblk, ONLY   : rn_efac 
    2121   USE sbcblk_phy     ! use some physical constants for flux computation 
    2222   ! 
     
    530530       
    531531      DO_2D_01_01 
    532          zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji  ,jj) + pssu(ji-1,jj) )   
    533          zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj  ) + pssv(ji,jj-1) )  
     532         zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * ( pssu(ji  ,jj) + pssu(ji-1,jj) )   
     533         zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * ( pssv(ji,jj  ) + pssv(ji,jj-1) )  
    534534      END_2D 
    535535      !  
     
    570570 
    571571#if defined key_si3 
    572          ! ------------------------------------------------------------ ! 
    573          !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
    574          ! ------------------------------------------------------------ ! 
    575          DO_2D_00_00 
    576              
    577             zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 
    578             zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 
    579     
    580             ptaui_ice(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             & 
    581                &                      +    rhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          & 
    582                &         * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 
    583             ptauj_ice(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1)             & 
    584                &                      +    rhoa(ji,jj  ) * pCd_du_ice(ji,jj  )  )          & 
    585                &         * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 
    586          END_2D 
    587          CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. ) 
    588          ! 
    589          IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   & 
    590             &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' ) 
     572      ! ------------------------------------------------------------ ! 
     573      !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
     574      ! ------------------------------------------------------------ ! 
     575      DO_2D_00_00            
     576         ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   & 
     577            &                      * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) ) 
     578         ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   & 
     579            &                      * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) ) 
     580      END_2D 
     581      CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. ) 
     582      ! 
     583      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   & 
     584         &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' ) 
    591585#endif 
    592586      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    599593   END SUBROUTINE abl_stp 
    600594!=================================================================================================== 
    601  
    602  
    603  
    604  
    605  
    606  
    607  
    608  
    609  
    610  
    611  
    612  
    613  
    614595 
    615596 
  • NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/src/ABL/sbcabl.F90

    r12551 r12565  
    335335      CALL blk_oce_1( kt,  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),   &   !   <<= in 
    336336         &                tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),   &   !   <<= in 
    337          &                sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m     ,   &   !   <<= in 
    338          &                sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) ,   &   !   <<= in 
     337         &                sf(jp_slp  )%fnow(:,:,1), sst_m, ssu_m, ssv_m     ,   &   !   <<= in 
     338         &                sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1),   &   !   <<= in 
     339         &                sf(jp_qsr  )%fnow(:,:,1), sf(jp_qlw  )%fnow(:,:,1),   &   !   <<= in 
    339340         &                tsk_m, zssq, zcd_du, zsen, zevp                       )   !   =>> out 
    340341 
  • NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/src/OCE/SBC/sbcblk.F90

    r12495 r12565  
    7474#endif 
    7575 
    76    INTEGER , PUBLIC            ::   jpfld         ! maximum number of files to read 
    77    INTEGER , PUBLIC, PARAMETER ::   jp_wndi = 1   ! index of 10m wind velocity (i-component) (m/s)    at T-point 
    78    INTEGER , PUBLIC, PARAMETER ::   jp_wndj = 2   ! index of 10m wind velocity (j-component) (m/s)    at T-point 
    79    INTEGER , PUBLIC, PARAMETER ::   jp_tair = 3   ! index of 10m air temperature             (Kelvin) 
    80    INTEGER , PUBLIC, PARAMETER ::   jp_humi = 4   ! index of specific humidity               ( % ) 
    81    INTEGER , PUBLIC, PARAMETER ::   jp_qsr  = 5   ! index of solar heat                      (W/m2) 
    82    INTEGER , PUBLIC, PARAMETER ::   jp_qlw  = 6   ! index of Long wave                       (W/m2) 
    83    INTEGER , PUBLIC, PARAMETER ::   jp_prec = 7   ! index of total precipitation (rain+snow) (Kg/m2/s) 
    84    INTEGER , PUBLIC, PARAMETER ::   jp_snow = 8   ! index of snow (solid prcipitation)       (kg/m2/s) 
    85    INTEGER , PUBLIC, PARAMETER ::   jp_slp  = 9   ! index of sea level pressure              (Pa) 
    86    INTEGER , PUBLIC, PARAMETER ::   jp_hpgi =10   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 
    87    INTEGER , PUBLIC, PARAMETER ::   jp_hpgj =11   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 
    88  
     76   INTEGER , PUBLIC, PARAMETER ::   jp_wndi  =  1   ! index of 10m wind velocity (i-component) (m/s)    at T-point 
     77   INTEGER , PUBLIC, PARAMETER ::   jp_wndj  =  2   ! index of 10m wind velocity (j-component) (m/s)    at T-point 
     78   INTEGER , PUBLIC, PARAMETER ::   jp_tair  =  3   ! index of 10m air temperature             (Kelvin) 
     79   INTEGER , PUBLIC, PARAMETER ::   jp_humi  =  4   ! index of specific humidity               ( % ) 
     80   INTEGER , PUBLIC, PARAMETER ::   jp_qsr   =  5   ! index of solar heat                      (W/m2) 
     81   INTEGER , PUBLIC, PARAMETER ::   jp_qlw   =  6   ! index of Long wave                       (W/m2) 
     82   INTEGER , PUBLIC, PARAMETER ::   jp_prec  =  7   ! index of total precipitation (rain+snow) (Kg/m2/s) 
     83   INTEGER , PUBLIC, PARAMETER ::   jp_snow  =  8   ! index of snow (solid prcipitation)       (kg/m2/s) 
     84   INTEGER , PUBLIC, PARAMETER ::   jp_slp   =  9   ! index of sea level pressure              (Pa) 
     85   INTEGER , PUBLIC, PARAMETER ::   jp_uoatm = 10   ! index of surface current (i-component) 
     86   !                                                !          seen by the atmospheric forcing (m/s) at T-point 
     87   INTEGER , PUBLIC, PARAMETER ::   jp_voatm = 11   ! index of surface current (j-component) 
     88   !                                                !          seen by the atmospheric forcing (m/s) at T-point 
     89   INTEGER , PUBLIC, PARAMETER ::   jp_hpgi  = 12   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 
     90   INTEGER , PUBLIC, PARAMETER ::   jp_hpgj  = 13   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 
     91   INTEGER , PUBLIC, PARAMETER ::   jpfld    = 13   ! maximum number of files to read 
     92 
     93   ! Warning: keep this structure allocatable for Agrif... 
    8994   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input atmospheric fields (file informations, fields read) 
    9095 
     
    98103   LOGICAL  ::   ln_Cd_L15      ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 
    99104   ! 
     105   LOGICAL  ::   ln_crt_fbk     ! Add surface current feedback to the wind stress computation  (Renault et al. 2020) 
     106   REAL(wp) ::   rn_stau_a      ! Alpha and Beta coefficients of Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta 
     107   REAL(wp) ::   rn_stau_b      !  
     108   ! 
    100109   REAL(wp)         ::   rn_pfac   ! multiplication factor for precipitation 
    101110   REAL(wp), PUBLIC ::   rn_efac   ! multiplication factor for evaporation 
    102    REAL(wp), PUBLIC ::   rn_vfac   ! multiplication factor for ice/ocean velocity in the calculation of wind stress 
    103111   REAL(wp)         ::   rn_zqt    ! z(q,t) : height of humidity and temperature measurements 
    104112   REAL(wp)         ::   rn_zu     ! z(u)   : height of wind measurements 
     
    162170      !! 
    163171      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files 
    164       TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i        ! array of namelist informations on the fields to read 
    165       TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read 
    166       TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        " 
    167       TYPE(FLD_N) ::   sn_slp , sn_hpgi, sn_hpgj               !       "                        " 
     172      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read 
     173      TYPE(FLD_N) ::   sn_wndi, sn_wndj , sn_humi, sn_qsr      ! informations about the fields to be read 
     174      TYPE(FLD_N) ::   sn_qlw , sn_tair , sn_prec, sn_snow     !       "                        " 
     175      TYPE(FLD_N) ::   sn_slp , sn_uoatm, sn_voatm             !       "                        " 
     176      TYPE(FLD_N) ::   sn_hpgi, sn_hpgj                        !       "                        " 
     177      INTEGER     ::   ipka                                    ! number of levels in the atmospheric variable 
    168178      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields 
    169          &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj,       & 
     179         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_uoatm, sn_voatm,     & 
     180         &                 sn_hpgi, sn_hpgj,                                          & 
    170181         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm 
    171182         &                 cn_dir , rn_zqt, rn_zu,                                    & 
    172          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           & 
     183         &                 rn_pfac, rn_efac, ln_Cd_L12, ln_Cd_L15,                    & 
     184         &                 ln_crt_fbk, rn_stau_a, rn_stau_b,                          &   ! current feedback 
    173185         &                 ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh  ! cool-skin / warm-layer !LB 
    174186      !!--------------------------------------------------------------------- 
     
    242254      !                                   !* set the bulk structure 
    243255      !                                      !- store namelist information in an array 
    244       IF( ln_blk ) jpfld = 9 
    245       IF( ln_abl ) jpfld = 11 
    246       ALLOCATE( slf_i(jpfld) ) 
    247       ! 
    248       slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj 
    249       slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw 
    250       slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi 
    251       slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow 
    252       slf_i(jp_slp ) = sn_slp 
    253       IF( ln_abl ) THEN 
    254          slf_i(jp_hpgi) = sn_hpgi   ;   slf_i(jp_hpgj) = sn_hpgj 
    255       END IF 
     256      ! 
     257      slf_i(jp_wndi ) = sn_wndi    ;   slf_i(jp_wndj ) = sn_wndj 
     258      slf_i(jp_qsr  ) = sn_qsr     ;   slf_i(jp_qlw  ) = sn_qlw 
     259      slf_i(jp_tair ) = sn_tair    ;   slf_i(jp_humi ) = sn_humi 
     260      slf_i(jp_prec ) = sn_prec    ;   slf_i(jp_snow ) = sn_snow 
     261      slf_i(jp_slp  ) = sn_slp 
     262      slf_i(jp_uoatm) = sn_uoatm   ;   slf_i(jp_voatm) = sn_voatm 
     263      slf_i(jp_hpgi ) = sn_hpgi    ;   slf_i(jp_hpgj ) = sn_hpgj 
     264      ! 
     265      IF( .NOT. ln_abl ) THEN   ! force to not use jp_hpgi and jp_hpgj, should already be done in namelist_* but we never know... 
     266         slf_i(jp_hpgi)%clname = 'NOT USED' 
     267         slf_i(jp_hpgj)%clname = 'NOT USED' 
     268      ENDIF 
    256269      ! 
    257270      !                                      !- allocate the bulk structure 
     
    264277      DO jfpr= 1, jpfld 
    265278         ! 
    266          IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to zero) 
    267             ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 
    268             sf(jfpr)%fnow(:,:,1) = 0._wp 
     279         IF(   ln_abl    .AND.                                                      & 
     280            &    ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR.   & 
     281            &      jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair     )  ) THEN 
     282            ipka = jpka   ! ABL: some fields are 3D input 
     283         ELSE 
     284            ipka = 1 
     285         ENDIF 
     286         ! 
     287         ALLOCATE( sf(jfpr)%fnow(jpi,jpj,ipka) ) 
     288         ! 
     289         IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to default) 
     290            IF(     jfpr == jp_slp  ) THEN 
     291               sf(jfpr)%fnow(:,:,1:ipka) = 101325._wp   ! use standard pressure in Pa 
     292            ELSEIF( jfpr == jp_prec .OR. jfpr == jp_snow .OR. jfpr == jp_uoatm .OR. jfpr == jp_voatm ) THEN 
     293               sf(jfpr)%fnow(:,:,1:ipka) = 0._wp        ! no precip or no snow or no surface currents 
     294            ELSEIF( ( jfpr == jp_hpgi .OR. jfpr == jp_hpgj ) .AND. .NOT. ln_abl ) THEN 
     295               DEALLOCATE( sf(jfpr)%fnow )              ! deallocate as not used in this case 
     296            ELSE 
     297               WRITE(ctmp1,*) 'sbc_blk_init: no default value defined for field number', jfpr 
     298               CALL ctl_stop( ctmp1 ) 
     299            ENDIF 
    269300         ELSE                                                  !-- used field  --! 
    270             IF(   ln_abl    .AND.                                                      & 
    271                &    ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR.   & 
    272                &      jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair     )  ) THEN   ! ABL: some fields are 3D input 
    273                ALLOCATE( sf(jfpr)%fnow(jpi,jpj,jpka) ) 
    274                IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,jpka,2) ) 
    275             ELSE                                                                                ! others or Bulk fields are 2D fiels 
    276                ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 
    277                IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) ) 
    278             ENDIF 
     301            IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,ipka,2) )   ! allocate array for temporal interpolation 
    279302            ! 
    280303            IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 )   & 
     
    327350         WRITE(numout,*) '      factor applied on precipitation (total & snow)      rn_pfac      = ', rn_pfac 
    328351         WRITE(numout,*) '      factor applied on evaporation                       rn_efac      = ', rn_efac 
    329          WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac 
    330352         WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))' 
    331353         WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12 
    332354         WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15 
     355         WRITE(numout,*) '      use surface current feedback on wind stress         ln_crt_fbk   = ', ln_crt_fbk 
     356         IF(ln_crt_fbk) THEN 
     357         WRITE(numout,*) '         Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta' 
     358         WRITE(numout,*) '            Alpha                                         rn_stau_a    = ', rn_stau_a 
     359         WRITE(numout,*) '            Beta                                          rn_stau_b    = ', rn_stau_b 
     360         ENDIF 
    333361         ! 
    334362         WRITE(numout,*) 
     
    429457      !                                            ! compute the surface ocean fluxes using bulk formulea 
    430458      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
    431          CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),   &   !   <<= in 
    432             &                sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),   &   !   <<= in 
    433             &                sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m,       &   !   <<= in 
    434             &                sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1),   &   !   <<= in (wl/cs) 
    435             &                tsk_m, zssq, zcd_du, zsen, zevp )                       !   =>> out 
    436  
    437          CALL blk_oce_2(     sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1),   &   !   <<= in 
    438             &                sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1),   &   !   <<= in 
    439             &                sf(jp_snow)%fnow(:,:,1), tsk_m,                     &   !   <<= in 
    440             &                zsen, zevp )                                            !   <=> in out 
     459         CALL blk_oce_1( kt, sf(jp_wndi )%fnow(:,:,1), sf(jp_wndj )%fnow(:,:,1),   &   !   <<= in 
     460            &                sf(jp_tair )%fnow(:,:,1), sf(jp_humi )%fnow(:,:,1),   &   !   <<= in 
     461            &                sf(jp_slp  )%fnow(:,:,1), sst_m, ssu_m, ssv_m,        &   !   <<= in 
     462            &                sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1),   &   !   <<= in 
     463            &                sf(jp_qsr  )%fnow(:,:,1), sf(jp_qlw  )%fnow(:,:,1),   &   !   <<= in (wl/cs) 
     464            &                tsk_m, zssq, zcd_du, zsen, zevp )                         !   =>> out 
     465 
     466         CALL blk_oce_2(     sf(jp_tair )%fnow(:,:,1), sf(jp_qsr  )%fnow(:,:,1),   &   !   <<= in 
     467            &                sf(jp_qlw  )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1),   &   !   <<= in 
     468            &                sf(jp_snow )%fnow(:,:,1), tsk_m,                      &   !   <<= in 
     469            &                zsen, zevp )                                              !   <=> in out 
    441470      ENDIF 
    442471      ! 
     
    470499 
    471500 
    472    SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, &  ! inp 
    473       &                  pslp , pst   , pu   , pv,        &  ! inp 
    474       &                  pqsr , pqlw  ,                   &  ! inp 
    475       &                  ptsk, pssq , pcd_du, psen , pevp   )  ! out 
     501   SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi,        &  ! inp 
     502      &                      pslp , pst  , pu   , pv,            &  ! inp 
     503      &                      puatm, pvatm, pqsr , pqlw ,         &  ! inp 
     504      &                      ptsk , pssq , pcd_du, psen, pevp   )   ! out 
    476505      !!--------------------------------------------------------------------- 
    477506      !!                     ***  ROUTINE blk_oce_1  *** 
     
    498527      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pu     ! surface current at U-point (i-component) [m/s] 
    499528      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pv     ! surface current at V-point (j-component) [m/s] 
     529      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   puatm  ! surface current seen by the atm at T-point (i-component) [m/s] 
     530      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pvatm  ! surface current seen by the atm at T-point (j-component) [m/s] 
    500531      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqsr   ! 
    501532      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqlw   ! 
     
    508539      INTEGER  ::   ji, jj               ! dummy loop indices 
    509540      REAL(wp) ::   zztmp                ! local variable 
     541      REAL(wp) ::   zstmax, zstau 
     542#if defined key_cyclone 
    510543      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
     544#endif 
     545      REAL(wp), DIMENSION(jpi,jpj) ::   ztau_i, ztau_j    ! wind stress components at T-point 
    511546      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    512547      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
     
    532567      zwnd_j(:,:) = 0._wp 
    533568      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012) 
    534       DO_2D_00_00 
    535          pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 
    536          pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 
     569      DO_2D_11_11 
     570         zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 
     571         zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 
     572         ! ... scalar wind at T-point (not masked) 
     573         wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) 
     574      END_2D 
     575#else 
     576      ! ... scalar wind module at T-point (not masked) 
     577      DO_2D_11_11 
     578         wndm(ji,jj) = SQRT(  pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj)  ) 
    537579      END_2D 
    538580#endif 
    539       DO_2D_00_00 
    540          zwnd_i(ji,jj) = (  pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
    541          zwnd_j(ji,jj) = (  pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
    542       END_2D 
    543       CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 
    544       ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
    545       wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   & 
    546          &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1) 
    547  
    548581      ! ----------------------------------------------------------------------------- ! 
    549582      !      I   Solar FLUX                                                           ! 
     
    597630         END_2D 
    598631      ENDIF 
    599  
    600  
    601632 
    602633      !! Time to call the user-selected bulk parameterization for 
     
    673704         pevp(:,:) = pevp(:,:) * tmask(:,:,1) 
    674705 
    675          ! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array... 
    676          zcd_oce = 0._wp 
    677          WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm 
    678          zwnd_i = zcd_oce * zwnd_i 
    679          zwnd_j = zcd_oce * zwnd_j 
    680  
    681          CALL iom_put( "taum_oce", taum )   ! output wind stress module 
     706         DO_2D_11_11 
     707            IF( wndm(ji,jj) > 0._wp ) THEN 
     708               zztmp = taum(ji,jj) / wndm(ji,jj) 
     709#if defined key_cyclone 
     710               ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
     711               ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
     712#else 
     713               ztau_i(ji,jj) = zztmp * pwndi(ji,jj) 
     714               ztau_j(ji,jj) = zztmp * pwndj(ji,jj) 
     715#endif 
     716            ELSE 
     717               ztau_i(ji,jj) = 0._wp 
     718               ztau_j(ji,jj) = 0._wp                  
     719            ENDIF 
     720         END_2D 
     721 
     722         IF( ln_crt_fbk ) THEN   ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715) 
     723            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) 
     724            DO_2D_01_01   ! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop 
     725               zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax )   ! stau (<0) must be smaller than zstmax 
     726               ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj  ) + pu(ji,jj) ) - puatm(ji,jj) ) 
     727               ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji  ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) ) 
     728               taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) ) 
     729            END_2D 
     730         ENDIF 
    682731 
    683732         ! ... utau, vtau at U- and V_points, resp. 
    684733         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
    685734         !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 
    686          DO_2D_10_10 
    687             utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) & 
    688                &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
    689             vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) & 
    690                &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
     735         DO_2D_00_00              ! start loop at 2, in case ln_crt_fbk = T 
     736            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj  ) ) & 
     737               &              * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
     738            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji  ,jj+1) ) & 
     739               &              * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
    691740         END_2D 
    692          CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
     741 
     742         IF( ln_crt_fbk ) THEN 
     743            CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1., taum, 'T', -1. ) 
     744         ELSE 
     745            CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
     746         ENDIF 
     747 
     748         CALL iom_put( "taum_oce", taum )   ! output wind stress module 
    693749 
    694750         IF(sn_cfctl%l_prtctl) THEN 
     
    861917      ! 
    862918      INTEGER  ::   ji, jj    ! dummy loop indices 
    863       REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point 
    864919      REAL(wp) ::   zootm_su                      ! sea-ice surface mean temperature 
    865920      REAL(wp) ::   zztmp1, zztmp2                ! temporary arrays 
     
    872927      ! ------------------------------------------------------------ ! 
    873928      ! C-grid ice dynamics :   U & V-points (same as ocean) 
    874       DO_2D_00_00 
    875          zwndi_t = (  pwndi(ji,jj) - rn_vfac * 0.5_wp * ( puice(ji-1,jj  ) + puice(ji,jj) )  ) 
    876          zwndj_t = (  pwndj(ji,jj) - rn_vfac * 0.5_wp * ( pvice(ji  ,jj-1) + pvice(ji,jj) )  ) 
    877          wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
     929      DO_2D_11_11 
     930         wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 
    878931      END_2D 
    879       CALL lbc_lnk( 'sbcblk', wndm_ice, 'T',  1. ) 
    880932      ! 
    881933      ! Make ice-atm. drag dependent on ice concentration 
     
    902954         ! C-grid ice dynamics :   U & V-points (same as ocean) 
    903955         DO_2D_00_00 
    904             putaui(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * zcd_dui(ji+1,jj)             & 
    905                &                      + rhoa(ji  ,jj) * zcd_dui(ji  ,jj)  )          & 
    906                &         * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) 
    907             pvtaui(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * zcd_dui(ji,jj+1)             & 
    908                &                      + rhoa(ji,jj  ) * zcd_dui(ji,jj  )  )          & 
    909                &         * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) 
     956            putaui(ji,jj) = 0.25_wp * ( rhoa(ji+1,jj) * zcd_dui(ji+1,jj) + rhoa(ji,jj) * zcd_dui(ji,jj) )   & 
     957               &                    * (                   pwndi(ji+1,jj) +                 pwndi(ji,jj) ) 
     958            pvtaui(ji,jj) = 0.25_wp * ( rhoa(ji,jj+1) * zcd_dui(ji,jj+1) + rhoa(ji,jj) * zcd_dui(ji,jj) )   & 
     959               &                    * (                   pwndj(ji,jj+1) +                 pwndj(ji,jj) ) 
    910960         END_2D 
    911961         CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. ) 
Note: See TracChangeset for help on using the changeset viewer.