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

Changeset 11286


Ignore:
Timestamp:
2019-07-18T11:54:22+02:00 (4 years ago)
Author:
dford
Message:

Merge in latest version of AMM15_v3_6_STABLE_package_collate.

Location:
branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM
Files:
29 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/CONFIG/SHARED/field_def.xml

    r5517 r11286  
    374374         <field id="uocet"        long_name="ocean transport along i-axis times temperature (CRS)"                                               unit="degC*m/s"   grid_ref="grid_U_3D" /> 
    375375         <field id="uoces"        long_name="ocean transport along i-axis times salinity (CRS)"                                                  unit="1e-3*m/s"   grid_ref="grid_U_3D" /> 
     376    <field id="ustokes"      long_name="Stokes Drift Velocity i-axis"                           standard_name="StokesDrift_x_velocity"      unit="m/s"        grid_ref="grid_U_3D" />   
     377    <field id="ustokes_e3u"  long_name="Stokes Drift Velocity i-axis  (thickness weighted)"                                                 unit="m/s"        grid_ref="grid_U_3D"  > ustokes * e3u </field> 
    376378 
    377379         <!-- variables available with MLE --> 
    … …  
    409411         <field id="vocet"        long_name="ocean transport along j-axis times temperature (CRS)"                                               unit="degC*m/s"   grid_ref="grid_V_3D" /> 
    410412         <field id="voces"        long_name="ocean transport along j-axis times salinity (CRS)"                                                  unit="1e-3*m/s"   grid_ref="grid_V_3D" /> 
     413    <field id="vstokes"      long_name="Stokes Drift Velocity j-axis"                           standard_name="StokesDrift_y_velocity"      unit="m/s"        grid_ref="grid_V_3D" />   
     414    <field id="vstokes_e3v"  long_name="Stokes Drift Velocity j-axis  (thickness weighted)"                                                 unit="m/s"        grid_ref="grid_V_3D"  > vstokes * e3v </field> 
    411415 
    412416         <!-- variables available with MLE --> 
    … …  
    437441        <field id="woce"         long_name="ocean vertical velocity"              standard_name="upward_sea_water_velocity"   unit="m/s"  /> 
    438442        <field id="wocetr_eff"   long_name="effective ocean vertical transport"                                               unit="m3/s" /> 
     443        <field id="wstokes"      long_name="Stokes Drift vertical velocity"       standard_name="upward_StokesDrift_velocity" unit="m/s" /> 
    439444 
    440445        <!-- woce_eiv: available with key_traldf_eiv and key_diaeiv --> 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/CONFIG/SHARED/namelist_ref

    r10728 r11286  
    240240   ln_blk_core = .true.    !  CORE bulk formulation                     (T => fill namsbc_core) 
    241241   ln_blk_mfs  = .false.   !  MFS bulk formulation                      (T => fill namsbc_mfs ) 
    242    ln_cpl      = .false.   !  atmosphere coupled   formulation          ( requires key_oasis3 ) 
    243    ln_mixcpl   = .false.   !  forced-coupled mixed formulation          ( requires key_oasis3 ) 
     242   ln_cpl      = .false.   !  coupled   formulation                       ( requires key_oasis3 )  
     243   ln_mixcpl   = .false.   !  forced-coupled mixed atmosphere formulation ( requires key_oasis3 )  
     244   ln_wavcpl   = .false.   !  forced-coupled mixed wave       formulation ( requires key_oasis3 ) 
    244245   nn_components = 0       !  configuration of the opa-sas OASIS coupling 
    245246                           !  =0 no opa-sas OASIS coupling: default single executable configuration 
    … …  
    264265                           !     =1 global mean of e-p-r set to zero at each time step 
    265266                           !     =2 annual global mean of e-p-r set to zero 
    266    ln_wave = .false.       !  Activate coupling with wave (either Stokes Drift or Drag coefficient, or both)  (T => fill namsbc_wave) 
    267    ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave) 
    268    ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave) 
     267   ln_wave = .false.       !  Activate coupling with wave (T => fill namsbc_wave)    
    269268   nn_lsm  = 0             !  =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) , 
    270269                           !  =1:n number of iterations of land/sea mask application for input fields (fill land/sea mask filename field) 
    … …  
    274273                           !  = 1  Average and redistribute per-category fluxes, forced mode only, not yet implemented coupled 
    275274                           !  = 2  Redistribute a single flux over categories (coupled mode only) 
     275   nn_drag   = 0      !  formula to calculate momentum from the wind components  
     276                           ! = 0 UKMO SHELF formulation  
     277                           ! = 1 standard formulation with forced of coupled drag coefficient  
     278                           ! = 2 standard formulation with constant drag coefficient  
     279                           ! = 3 momentum calculated from core forcing fields 
    276280/ 
    277281!----------------------------------------------------------------------- 
    … …  
    296300   sn_emp      = 'emp'       ,        24         , 'emp'     , .false.      , .false., 'yearly'  , ''       , ''       , '' 
    297301 
    298    cn_dir      = './'      !  root directory for the location of the flux files 
     302   cn_dir       = './'      !  root directory for the location of the flux files  
     303   ln_shelf_flx = .false.   !  UKMO SHELF specific flux flag - read from file wind components instead of momentum   
     304   ln_rel_wind  = .false.   !  UKMO SHELF - calculate momentum from wind speed relative to currents   
     305   rn_wfac      = 1.0       !  UKMO SHELF - multiplicative factor for ocean/wind velocity 
    299306/ 
    300307!----------------------------------------------------------------------- 
    … …  
    363370   sn_snd_crt    =       'none'                 ,    'no'    , 'spherical' , 'eastward-northward' ,  'T' 
    364371   sn_snd_co2    =       'coupled'              ,    'no'    ,     ''      ,         ''           ,   '' 
     372   sn_snd_crtw   =       'none'                 ,    'no'    ,     ''      ,         ''           , 'U,V'   
     373   sn_snd_ifrac  =       'none'                 ,    'no'    ,     ''      ,         ''           ,   ''   
     374   sn_snd_wlev   =       'none'                 ,    'no'    ,     ''      ,         ''           ,   '' 
    365375! receive 
    366376   sn_rcv_w10m   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
    … …  
    374384   sn_rcv_cal    =       'coupled'              ,    'no'    ,     ''      ,         ''          ,   '' 
    375385   sn_rcv_co2    =       'coupled'              ,    'no'    ,     ''      ,         ''          ,   '' 
     386   sn_rcv_hsig   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
     387   sn_rcv_iceflx =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
     388   sn_rcv_mslp   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
     389   sn_rcv_phioc  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
     390   sn_rcv_sdrft  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
     391   sn_rcv_wper   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
     392   sn_rcv_wfreq  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
     393   sn_rcv_wnum   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
     394   sn_rcv_tauoc  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
     395   sn_rcv_tauw   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
     396   sn_rcv_wdrag  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
    376397! 
    377398   nn_cplmodel   =     1     !  Maximum number of models to/from which NEMO is potentialy sending/receiving data 
    378399   ln_usecplmask = .false.   !  use a coupling mask file to merge data received from several models 
    379400                             !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel) 
     401   ln_usernfmask = .false.   !  use a runoff mask file to merge data received from several models  
     402                             !   -> file rnfmask.nc with the float variable called rnfmask (jpi,jpj,nn_cplmodel) 
    380403/ 
    381404!----------------------------------------------------------------------- 
    … …  
    9841007   rn_clim_galp  = 0.267   !  galperin limit 
    9851008   ln_sigpsi     = .true.  !  Activate or not Burchard 2001 mods on psi schmidt number in the wb case 
    986    rn_crban      = 100.    !  Craig and Banner 1994 constant for wb tke flux 
     1009   rn_crban_default = 100.    !  Craig and Banner 1994 constant for wb tke flux 
    9871010   rn_charn      = 70000.  !  Charnock constant for wb induced roughness length 
    9881011   rn_hsro       =  0.02   !  Minimum surface roughness 
    9891012   rn_frac_hs    =   1.3   !  Fraction of wave height as roughness (if nn_z0_met=2) 
    990    nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2) 
     1013   nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2/3)  
     1014      !                             ! =3 requires ln_wave=T  
    9911015   nn_bc_surf    =     1   !  surface condition (0/1=Dir/Neum) 
    9921016   nn_bc_bot     =     1   !  bottom condition (0/1=Dir/Neum) 
    … …  
    13001324/ 
    13011325!----------------------------------------------------------------------- 
    1302 &namsbc_wave   ! External fields from wave model 
     1326&namsbc_wave   ! External fields from wave model                        (ln_wave=T) 
    13031327!----------------------------------------------------------------------- 
    13041328!              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
    13051329!              !             !  (if <0  months)  !   name       !   (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      ! 
    1306    sn_cdg      =  'cdg_wave' ,        1          , 'drag_coeff' ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
     1330   ln_cdgw     = .false.   !  Neutral drag coefficient read from wave model  
     1331   ln_sdw      = .false.   !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift  
     1332   ln_stcor    = .false.   !  Activate Stokes Coriolis term  
     1333   ln_tauoc    = .false.   !  Activate ocean stress modified by external wave induced stress  
     1334   ln_tauw     = .false.   !  Activate ocean stress components from wave model  
     1335   ln_phioc    = .false.   !  Activate wave to ocean energy  
     1336   ln_rough    = .false.   !  Wave roughness equals the significant wave height  
     1337   ln_zdfqiao  = .false.   !  Enhanced wave vertical mixing Qiao (2010)  
     1338   nn_sdrift   =  0   !  Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift  
     1339                           ! = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]  
     1340                           ! = 1 Phillips:                      v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))]  
     1341   nn_wmix     =  0   !  type of wave breaking mixing  
     1342                           ! = 0 Craig and Banner formulation (original NEMO formulation)   
     1343                           ! = 1 Janssen formulation (no assumption of direct energy conversion) 
     1344   sn_cdg      =  'cdw_wave' ,        1          , 'drag_coeff' ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    13071345   sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    13081346   sn_vsd      =  'sdw_wave' ,        1          , 'v_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    1309    sn_wn       =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    1310 ! 
    1311    cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files 
     1347   sn_hsw      =  'sdw_wave' ,        1          , 'hs'         ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
     1348   sn_wmp      =  'sdw_wave' ,        1          , 'wmp'        ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
     1349   sn_wfr      =  'sdw_wave' ,        1          , 'wave_freq'  ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
     1350   sn_wnum     =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
     1351   sn_tauoc    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
     1352   sn_tauwx    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
     1353   sn_tauwy    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
     1354   sn_phioc    =  'sdw_wave' ,        1          , 'wave_energy',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
     1355   cn_dir      = './'  !  root directory for the location of wave forcing files 
     1356!   
    13121357/ 
    13131358!----------------------------------------------------------------------- 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/CONFIG/cfg.txt

    r8058 r11286  
    1111ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
    1212ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 
     13GYRE_LONG OPA_SRC 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/DIA/diadimg.F90

    r8058 r11286  
    124124 
    125125    CASE DEFAULT 
    126        IF(lwp) WRITE(numout,*) ' E R R O R : bad cd_type in dia_wri_dimg ' 
    127        STOP 'dia_wri_dimg' 
    128  
     126       WRITE(numout,*) ' E R R O R : bad cd_type in dia_wri_dimg'  
     127       CALL ctl_stop( 'STOP', 'dia_wri_dimg :bad cd_type in dia_wri_dimg ' )  
    129128    END SELECT 
    130129 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r8561 r11286  
    4141   USE zdfddm          ! vertical  physics: double diffusion 
    4242   USE diahth          ! thermocline diagnostics 
     43   USE sbcwave         ! wave parameters 
    4344   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    4445   USE in_out_manager  ! I/O manager 
    … …  
    709710         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un 
    710711            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout ) 
     712         IF( ln_wave .AND. ln_sdw) THEN   
     713            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd   
     714               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )   
     715         ENDIF  
    711716         IF( ln_traldf_gdia ) THEN 
    712717            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv 
    … …  
    727732         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn 
    728733            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout ) 
     734         IF( ln_wave .AND. ln_sdw) THEN   
     735            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd   
     736               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )   
     737         ENDIF   
    729738         IF( ln_traldf_gdia ) THEN 
    730739            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv 
    … …  
    763772               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 
    764773         ENDIF 
     774            
     775         IF( ln_wave .AND. ln_sdw) THEN   
     776            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd   
     777               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )   
     778         ENDIF   
    765779         !                                                                                      !!! nid_W : 2D 
    766780#if defined key_traldf_c2d 
    … …  
    943957#endif 
    944958 
     959      IF( ln_wave .AND. ln_sdw ) THEN   
     960         CALL histwrite( nid_U, "sdzocrtx", it, usd           , ndim_U , ndex_U )    ! i-StokesDrift-current   
     961         CALL histwrite( nid_V, "sdmecrty", it, vsd           , ndim_V , ndex_V )    ! j-StokesDrift-current   
     962         CALL histwrite( nid_W, "sdvecrtz", it, wsd           , ndim_T , ndex_T )    ! StokesDrift vert. current   
     963      ENDIF   
     964        
    945965      ! 3. Close all files 
    946966      ! --------------------------------------- 
    … …  
    10481068            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
    10491069      END IF 
     1070      !   
     1071      IF( ln_wave .AND. ln_sdw ) THEN   
     1072         CALL histdef( id_i, "sdzocrtx", "Stokes Drift Zonal", "m/s"    , &   ! StokesDrift zonal current   
     1073            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )   
     1074         CALL histdef( id_i, "sdmecrty", "Stokes Drift Merid", "m/s"    , &   ! StokesDrift meridonal current   
     1075            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )   
     1076         CALL histdef( id_i, "sdvecrtz", "Stokes Drift Vert", "m/s"    , &   ! StokesDrift vertical current   
     1077            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )   
     1078      ENDIF   
    10501079 
    10511080#if defined key_lim2 
    … …  
    10651094 
    10661095      ! Write all fields on T grid 
     1096      IF( ln_wave .AND. ln_sdw ) THEN   
     1097         CALL histwrite( id_i, "sdzocrtx", kt, usd           , jpi*jpj*jpk, idex)     ! now StokesDrift i-velocity   
     1098         CALL histwrite( id_i, "sdmecrty", kt, vsd           , jpi*jpj*jpk, idex)     ! now StokesDrift j-velocity   
     1099         CALL histwrite( id_i, "sdvecrtz", kt, wsd           , jpi*jpj*jpk, idex)     ! now StokesDrift k-velocity   
     1100      ENDIF 
    10671101      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature 
    10681102      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r10268 r11286  
    1111   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping 
    1212   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility 
     13   !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence 
    1314   !!--------------------------------------------------------------------- 
    1415#if defined key_dynspg_ts   ||   defined key_esopa 
    … …  
    3132   USE sbctide         ! tides 
    3233   USE updtide         ! tide potential 
     34   USE sbcwave         ! surface wave  
     35   !  
    3336   USE lib_mpp         ! distributed memory computing library 
    3437   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    … …  
    333336         END DO 
    334337      END DO 
     338               
     339      !!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....   
     340      !!gm  Is it correct to do so ?   I think so...   
     341               
     342               
    335343      !                                   !* barotropic Coriolis trends (vorticity scheme dependent) 
    336344      !                                   ! -------------------------------------------------------- 
    … …  
    475483      ! 
    476484      ! ----------------------------------------------------------------------- 
    477       !  Phase 2 : Integration of the barotropic equations  
     485      !  Phase 2 : Integration of the barotropic equations 
    478486      ! ----------------------------------------------------------------------- 
    479487      ! 
    480488      !                                             ! ==================== ! 
    481489      !                                             !    Initialisations   ! 
    482       !                                             ! ==================== !   
     490      !                                             ! ==================== ! 
    483491      ! Initialize barotropic variables:       
    484492      IF( ll_init )THEN 
    … …  
    489497         ub_e   (:,:) = 0._wp 
    490498         vb_e   (:,:) = 0._wp 
     499      ENDIF 
     500      !  
     501      ! Moved here to allow merging with other branch   
     502      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary   
     503         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)   
    491504      ENDIF 
    492505      ! 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r8058 r11286  
    1616   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
    1717   !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity  
     18   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T) 
    1819   !!---------------------------------------------------------------------- 
    1920 
    … …  
    3233   USE trd_oce        ! trends: ocean variables 
    3334   USE trddyn         ! trend manager: dynamics 
     35   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)  
     36   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force  
     37   !  
    3438   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3539   USE prtctl         ! Print control 
    … …  
    9195      ! 
    9296      CASE ( -1 )                                      ! esopa: test all possibility with control print 
    93          CALL vor_ene( kt, ntot, ua, va ) 
     97         CALL vor_ene( kt, ntot, un, vn, ua, va ) 
    9498         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, & 
    9599            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    96          CALL vor_ens( kt, ntot, ua, va ) 
     100         CALL vor_ens( kt, ntot, un, vn, ua, va ) 
    97101         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, & 
    98102            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    … …  
    100104         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, & 
    101105            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    102          CALL vor_een( kt, ntot, ua, va ) 
     106         CALL vor_een( kt, ntot, un, vn, ua, va ) 
    103107         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, & 
    104108            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    … …  
    108112            ztrdu(:,:,:) = ua(:,:,:) 
    109113            ztrdv(:,:,:) = va(:,:,:) 
    110             CALL vor_ene( kt, nrvm, ua, va )                ! relative vorticity or metric trend 
     114            CALL vor_ene( kt, nrvm, un, vn, ua, va )        ! relative vorticity or metric trend 
    111115            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    112116            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    … …  
    114118            ztrdu(:,:,:) = ua(:,:,:) 
    115119            ztrdv(:,:,:) = va(:,:,:) 
    116             CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend 
     120            CALL vor_ene( kt, ncor, un, vn, ua, va )        ! planetary vorticity trend 
    117121            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    118122            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    119123            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    120124         ELSE 
    121             CALL vor_ene( kt, ntot, ua, va )                ! total vorticity 
     125                             CALL vor_ene( kt, ntot, un, vn, ua, va )    ! total vorticity trend  
     126            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend  
    122127         ENDIF 
    123128         ! 
    … …  
    126131            ztrdu(:,:,:) = ua(:,:,:) 
    127132            ztrdv(:,:,:) = va(:,:,:) 
    128             CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend 
     133            CALL vor_ens( kt, nrvm, un, vn, ua, va )        ! relative vorticity or metric trend 
    129134            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    130135            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    … …  
    132137            ztrdu(:,:,:) = ua(:,:,:) 
    133138            ztrdv(:,:,:) = va(:,:,:) 
    134             CALL vor_ens( kt, ncor, ua, va )                ! planetary vorticity trend 
     139            CALL vor_ens( kt, ncor, un, vn, ua, va )        ! planetary vorticity trend 
    135140            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    136141            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    137142            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    138143         ELSE 
    139             CALL vor_ens( kt, ntot, ua, va )                ! total vorticity 
     144                             CALL vor_ens( kt, ntot, un, vn, ua, va )    ! total vorticity  
     145            IF( ln_stcor )   CALL vor_ens( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend 
    140146         ENDIF 
    141147         ! 
    … …  
    144150            ztrdu(:,:,:) = ua(:,:,:) 
    145151            ztrdv(:,:,:) = va(:,:,:) 
    146             CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens) 
     152            CALL vor_ens( kt, nrvm, un, vn, ua, va )        ! relative vorticity or metric trend (ens) 
    147153            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    148154            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    … …  
    150156            ztrdu(:,:,:) = ua(:,:,:) 
    151157            ztrdv(:,:,:) = va(:,:,:) 
    152             CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene) 
     158            CALL vor_ene( kt, ncor, un, vn, ua, va )        ! planetary vorticity trend (ene) 
    153159            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    154160            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    155161            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    156162         ELSE 
    157             CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene) 
     163                             CALL vor_ens( kt, nrvm, un , vn , ua, va )   ! relative vorticity or metric trend (ens)  
     164                             CALL vor_ene( kt, ncor, un , vn , ua, va )   ! planetary vorticity trend (ene)  
     165            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    158166         ENDIF 
    159167         ! 
    … …  
    162170            ztrdu(:,:,:) = ua(:,:,:) 
    163171            ztrdv(:,:,:) = va(:,:,:) 
    164             CALL vor_een( kt, nrvm, ua, va )                ! relative vorticity or metric trend 
     172            CALL vor_een( kt, nrvm, un, vn, ua, va )        ! relative vorticity or metric trend 
    165173            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    166174            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    … …  
    168176            ztrdu(:,:,:) = ua(:,:,:) 
    169177            ztrdv(:,:,:) = va(:,:,:) 
    170             CALL vor_een( kt, ncor, ua, va )                ! planetary vorticity trend 
     178            CALL vor_een( kt, ncor, un, vn, ua, va )        ! planetary vorticity trend 
    171179            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    172180            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    173181            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    174182         ELSE 
    175             CALL vor_een( kt, ntot, ua, va )                ! total vorticity 
     183                             CALL vor_een( kt, ntot, un, vn, ua, va )    ! total vorticity  
     184            IF( ln_stcor )   CALL vor_een( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend 
    176185         ENDIF 
    177186         ! 
    … …  
    189198 
    190199 
    191    SUBROUTINE vor_ene( kt, kvor, pua, pva ) 
     200   SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva ) 
    192201      !!---------------------------------------------------------------------- 
    193202      !!                  ***  ROUTINE vor_ene  *** 
    … …  
    214223      !!---------------------------------------------------------------------- 
    215224      ! 
    216       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    217       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
    218       !                                                           ! =nrvm (relative vorticity or metric) 
    219       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    220       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
     225      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index  
     226      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;  
     227      !                                                                ! =nrvm (relative vorticity or metric)  
     228      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua         ! total u-trend  
     229      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva         ! total v-trend  
     230      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities 
    221231      ! 
    222232      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    … …  
    250260            DO jj = 1, jpjm1 
    251261               DO ji = 1, fs_jpim1   ! vector opt. 
    252                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    253                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     262                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &  
     263                       &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    254264                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) 
    255265               END DO 
    … …  
    260270               DO ji = 1, fs_jpim1   ! vector opt. 
    261271                  zwz(ji,jj) = ( ff (ji,jj)                                                                       & 
    262                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    263                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     272                       &       + (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &  
     273                       &           - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &  
    264274                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               & 
    265275                       &       ) 
    … …  
    270280         IF( ln_sco ) THEN 
    271281            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk) 
    272             zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 
    273             zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
     282            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * pun(:,:,jk)  
     283            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * pvn(:,:,jk) 
    274284         ELSE 
    275             zwx(:,:) = e2u(:,:) * un(:,:,jk) 
    276             zwy(:,:) = e1v(:,:) * vn(:,:,jk) 
     285            zwx(:,:) = e2u(:,:) * pun(:,:,jk) 
     286            zwy(:,:) = e1v(:,:) * pvn(:,:,jk) 
    277287         ENDIF 
    278288 
    … …  
    418428 
    419429 
    420    SUBROUTINE vor_ens( kt, kvor, pua, pva ) 
     430   SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva ) 
    421431      !!---------------------------------------------------------------------- 
    422432      !!                ***  ROUTINE vor_ens  *** 
    … …  
    443453      !!---------------------------------------------------------------------- 
    444454      ! 
    445       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    446       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
    447          !                                                        ! =nrvm (relative vorticity or metric) 
    448       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    449       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
     455      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index  
     456      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;  
     457         !                                                             ! =nrvm (relative vorticity or metric)  
     458      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua         ! total u-trend  
     459      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva         ! total v-trend  
     460      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities 
    450461      ! 
    451462      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    … …  
    479490            DO jj = 1, jpjm1 
    480491               DO ji = 1, fs_jpim1   ! vector opt. 
    481                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    482                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     492                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &  
     493                       &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    483494                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) 
    484495               END DO 
    … …  
    489500               DO ji = 1, fs_jpim1   ! vector opt. 
    490501                  zwz(ji,jj) = ( ff (ji,jj)                                                                       & 
    491                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    492                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     502                       &       + (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &  
     503                       &           - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    493504                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                & 
    494505                       &       ) 
    … …  
    501512               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking 
    502513                  zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk) 
    503                   zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) 
    504                   zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) 
     514                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 
     515                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 
    505516               END DO 
    506517            END DO 
    … …  
    508519            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop  
    509520               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking 
    510                   zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk) 
    511                   zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk) 
     521                  zwx(ji,jj) = e2u(ji,jj) * pun(ji,jj,jk) 
     522                  zwy(ji,jj) = e1v(ji,jj) * pvn(ji,jj,jk) 
    512523               END DO 
    513524            END DO 
    … …  
    536547 
    537548 
    538    SUBROUTINE vor_een( kt, kvor, pua, pva ) 
     549   SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva ) 
    539550      !!---------------------------------------------------------------------- 
    540551      !!                ***  ROUTINE vor_een  *** 
    … …  
    554565      !!---------------------------------------------------------------------- 
    555566      ! 
    556       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    557       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
    558       !                                                           ! =nrvm (relative vorticity or metric) 
    559       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    560       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
     567      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index  
     568      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ;  
     569      !                                                                ! =nrvm (relative vorticity or metric)  
     570      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua         ! total u-trend  
     571      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva         ! total v-trend  
     572      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities 
    561573      !! 
    562574      INTEGER  ::   ji, jj, jk                                    ! dummy loop indices 
    … …  
    604616                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
    605617                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) 
    606                      IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = 4.0_wp / ze3 
     618                     IF( ze3 /= 0._wp ) THEN  ;   ze3f(ji,jj,jk) = 4.0_wp / ze3  
     619                     ELSE                     ;   ze3f(ji,jj,jk) = 0._wp  
     620                     ENDIF 
    607621                  END DO 
    608622               END DO 
    … …  
    616630                     zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   & 
    617631                        &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) ) 
    618                      IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3 
     632                     IF( ze3 /= 0._wp ) THEN  ;  ze3f(ji,jj,jk) = zmsk / ze3  
     633                     ELSE                     ;  ze3f(ji,jj,jk) = 0._wp  
     634                     ENDIF 
    619635                  END DO 
    620636               END DO 
    … …  
    643659            DO jj = 1, jpjm1 
    644660               DO ji = 1, fs_jpim1   ! vector opt. 
    645                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    646                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     661                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &  
     662                       &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    647663                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk) 
    648664               END DO 
    … …  
    655671               DO ji = 1, fs_jpim1   ! vector opt. 
    656672                  zwz(ji,jj) = ( ff (ji,jj)                                                                       & 
    657                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    658                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     673                       &       + (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &  
     674                       &           - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    659675                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                & 
    660676                       &       ) * ze3f(ji,jj,jk) 
    … …  
    664680         END SELECT 
    665681 
    666          zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 
    667          zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
     682         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * pun(:,:,jk) 
     683         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * pvn(:,:,jk) 
    668684 
    669685         ! Compute and add the vorticity term trend 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90

    r8058 r11286  
    6262   USE lbcnfd         ! north fold treatment 
    6363   USE in_out_manager ! I/O manager 
     64   USE mod_oasis      ! coupling routines 
    6465 
    6566   IMPLICIT NONE 
    … …  
    20232024      !!---------------------------------------------------------------------- 
    20242025      ! 
     2026#if defined key_oasis3 || defined key_oasis3mct    
     2027      ! If we're trying to shut down cleanly then we need to consider the fact    
     2028      ! that this could be part of an MPMD configuration - we don't want to    
     2029      ! leave other components deadlocked.    
     2030      CALL oasis_abort(nproc,"mppstop","NEMO initiated abort")    
     2031#else    
    20252032      CALL mppsync 
    20262033      CALL mpi_finalize( info ) 
     2034#endif  
    20272035      ! 
    20282036   END SUBROUTINE mppstop 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90

    r8058 r11286  
    2424   !!   cpl_finalize : finalize the coupled mode communication 
    2525   !!---------------------------------------------------------------------- 
    26 #if defined key_oasis3 
     26#if defined key_oasis3 || defined key_oasis3mct 
    2727   USE mod_oasis                    ! OASIS3-MCT module 
    2828#endif 
    … …  
    3232   USE lbclnk                       ! ocean lateral boundary conditions (or mpp link) 
    3333 
     34#if defined key_cpl_rootexchg    
     35   USE lib_mpp, only : mppsync    
     36   USE lib_mpp, only : mppscatter,mppgather    
     37#endif     
     38 
    3439   IMPLICIT NONE 
    3540   PRIVATE 
    … …  
    4146   PUBLIC   cpl_freq 
    4247   PUBLIC   cpl_finalize 
     48#if defined key_mpp_mpi    
     49   INCLUDE 'mpif.h'    
     50#endif    
     51       
     52   INTEGER, PARAMETER         :: localRoot  = 0    
     53   LOGICAL                    :: commRank            ! true for ranks doing OASIS communication    
     54#if defined key_cpl_rootexchg    
     55   LOGICAL                    :: rootexchg =.true.   ! logical switch     
     56#else    
     57   LOGICAL                    :: rootexchg =.false.  ! logical switch     
     58#endif     
    4359 
    4460   INTEGER, PUBLIC            ::   OASIS_Rcv  = 1    !: return code if received field 
    … …  
    4662   INTEGER                    ::   ncomp_id          ! id returned by oasis_init_comp 
    4763   INTEGER                    ::   nerror            ! return error code 
    48 #if ! defined key_oasis3 
     64#if ! defined key_oasis3 && ! defined key_oasis3mct 
    4965   ! OASIS Variables not used. defined only for compilation purpose 
    5066   INTEGER                    ::   OASIS_Out         = -1 
    … …  
    6581   INTEGER                    ::   nsnd         ! total number of fields sent  
    6682   INTEGER                    ::   ncplmodel    ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 
    67    INTEGER, PUBLIC, PARAMETER ::   nmaxfld=50   ! Maximum number of coupling fields 
     83   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=55   ! Maximum number of coupling fields 
    6884   INTEGER, PUBLIC, PARAMETER ::   nmaxcat=5    ! Maximum number of coupling fields 
    6985   INTEGER, PUBLIC, PARAMETER ::   nmaxcpl=5    ! Maximum number of coupling fields 
    … …  
    8298 
    8399   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   exfld   ! Temporary buffer for receiving 
     100   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   tbuf  ! Temporary buffer for sending / receiving     
     101   INTEGER, PUBLIC :: localComm     
    84102 
    85103   !!---------------------------------------------------------------------- 
    … …  
    120138      IF ( nerror /= OASIS_Ok ) & 
    121139         CALL oasis_abort (ncomp_id, 'cpl_init','Failure in oasis_get_localcomm' ) 
     140      localComm = kl_comm 
    122141      ! 
    123142   END SUBROUTINE cpl_init 
    … …  
    149168      IF(lwp) WRITE(numout,*) 
    150169 
     170      commRank = .false.    
     171      IF ( rootexchg ) THEN    
     172         IF ( nproc == localRoot ) commRank = .true.    
     173      ELSE    
     174         commRank = .true.    
     175      ENDIF    
     176 
    151177      ncplmodel = kcplmodel 
    152178      IF( kcplmodel > nmaxcpl ) THEN 
    … …  
    172198      ishape(:,2) = (/ 1, nlej-nldj+1 /) 
    173199      ! 
    174       ! ... Allocate memory for data exchange 
    175       ! 
    176       ALLOCATE(exfld(nlei-nldi+1, nlej-nldj+1), stat = nerror) 
    177       IF( nerror > 0 ) THEN 
    178          CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating exfld')   ;   RETURN 
    179       ENDIF 
    180200      ! 
    181201      ! ----------------------------------------------------------------- 
    … …  
    183203      ! ----------------------------------------------------------------- 
    184204       
     205      IF ( rootexchg ) THEN    
     206        
     207         paral(1) = 2              ! box partitioning    
     208         paral(2) = 0              ! NEMO lower left corner global offset         
     209         paral(3) = jpiglo         ! local extent in i     
     210         paral(4) = jpjglo         ! local extent in j    
     211         paral(5) = jpiglo         ! global extent in x    
     212        
     213      ELSE     
    185214      paral(1) = 2                                              ! box partitioning 
    186215      paral(2) = jpiglo * (nldj-1+njmpp-1) + (nldi-1+nimpp-1)   ! NEMO lower left corner global offset     
    … …  
    196225      ENDIF 
    197226       
    198       CALL oasis_def_partition ( id_part, paral, nerror ) 
     227      ENDIF  
     228      IF ( commRank )  CALL oasis_def_partition ( id_part, paral, nerror )    
     229        
     230      ! ... Allocate memory for data exchange    
     231      !    
     232      ALLOCATE(exfld(paral(3), paral(4)), stat = nerror)    
     233      IF( nerror > 0 ) THEN    
     234         CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating exfld')   ;   RETURN    
     235      ENDIF    
     236      IF ( rootexchg ) THEN    
     237         ! Should possibly use one of the work arrays for tbuf really    
     238         ALLOCATE(tbuf(jpi, jpj, jpnij), stat = nerror)    
     239         IF( nerror > 0 ) THEN    
     240            CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating tbuf') ; RETURN    
     241         ENDIF    
     242      ENDIF                 
     243      !    
     244      IF (commRank ) THEN  
    199245      ! 
    200246      ! ... Announce send variables.  
    … …  
    288334         ENDIF 
    289335      END DO 
     336      !   
     337      ENDIF  ! commRank=true   
    290338       
    291339      !------------------------------------------------------------------ 
    … …  
    293341      !------------------------------------------------------------------ 
    294342       
    295       CALL oasis_enddef(nerror) 
    296       IF( nerror /= OASIS_Ok )   CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in oasis_enddef') 
     343      IF ( commRank ) THEN         
     344         CALL oasis_enddef(nerror)    
     345         IF( nerror /= OASIS_Ok )   CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in oasis_enddef')    
     346      ENDIF    
    297347      ! 
    298348   END SUBROUTINE cpl_define 
    … …  
    311361      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdata 
    312362      !! 
    313       INTEGER                                   ::   jc,jm     ! local loop index 
     363      INTEGER                                   ::   jn,jc,jm  ! local loop index 
    314364      !!-------------------------------------------------------------------- 
    315365      ! 
    … …  
    320370         
    321371            IF( ssnd(kid)%nid(jc,jm) /= -1 ) THEN 
    322                CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(nldi:nlei, nldj:nlej,jc), kinfo ) 
     372               IF ( rootexchg ) THEN    
     373                  !    
     374                  ! collect data on the local root process    
     375                  !    
     376                  CALL mppgather (pdata(:,:,jc),localRoot,tbuf)     
     377                  CALL mppsync     
     378                          
     379                  IF ( nproc == localRoot ) THEN    
     380                     DO jn = 1, jpnij    
     381                        exfld(nimppt(jn)-1+nldit(jn):nimppt(jn)+nleit(jn)-1,njmppt(jn)-1+nldjt(jn):njmppt(jn)+nlejt(jn)-1)= &    
     382                          tbuf(nldit(jn):nleit(jn),nldjt(jn):nlejt(jn),jn)    
     383                     ENDDO    
     384                     ! snd data to OASIS3    
     385                     CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, exfld, kinfo )    
     386                  ENDIF    
     387               ELSE    
     388                  ! snd data to OASIS3    
     389                  CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(nldi:nlei, nldj:nlej,jc), kinfo )    
     390               ENDIF    
    323391                
    324392               IF ( ln_ctl ) THEN         
    … …  
    358426      INTEGER                   , INTENT(  out) ::   kinfo     ! OASIS3 info argument 
    359427      !! 
    360       INTEGER                                   ::   jc,jm     ! local loop index 
     428      INTEGER                                   ::   jn,jc,jm  ! local loop index 
    361429      LOGICAL                                   ::   llaction, llfisrt 
    362430      !!-------------------------------------------------------------------- 
    … …  
    372440 
    373441            IF( srcv(kid)%nid(jc,jm) /= -1 ) THEN 
    374  
    375                CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, exfld, kinfo )          
     442               !    
     443               ! receive data from OASIS3    
     444               !    
     445               IF ( commRank )   CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, exfld, kinfo )    
     446               IF ( rootexchg )  CALL MPI_BCAST ( kinfo, 1, MPI_INTEGER, localRoot, localComm, nerror )  
    376447                
    377448               llaction =  kinfo == OASIS_Recvd   .OR. kinfo == OASIS_FromRest .OR.   & 
    … …  
    384455                  kinfo = OASIS_Rcv 
    385456                  IF( llfisrt ) THEN  
    386                      pdata(nldi:nlei,nldj:nlej,jc) =                                 exfld(:,:) * pmask(nldi:nlei,nldj:nlej,jm) 
     457                     IF ( rootexchg ) THEN    
     458                        ! distribute data to processes    
     459                        !    
     460                        IF ( nproc == localRoot ) THEN    
     461                           DO jn = 1, jpnij    
     462                              tbuf(nldit(jn):nleit(jn),nldjt(jn):nlejt(jn),jn)=          &    
     463                              exfld(nimppt(jn)-1+nldit(jn):nimppt(jn)+nleit(jn)-1,njmppt(jn)-1+nldjt(jn):njmppt(jn)+nlejt(jn)-1)    
     464                              ! NOTE: we are missing combining this with pmask (see else below)    
     465                           ENDDO    
     466                        ENDIF    
     467                        CALL mppscatter(tbuf,localRoot,pdata(:,:,jc))     
     468                        CALL mppsync    
     469                     ELSE    
     470                        pdata(nldi:nlei, nldj:nlej, jc) = exfld(:,:) * pmask(nldi:nlei,nldj:nlej,jm)    
     471                     ENDIF    
    387472                     llfisrt = .FALSE. 
    388473                  ELSE 
    … …  
    462547#if defined key_oa3mct_v3 
    463548         CALL oasis_get_freqs(id, mop, 1, itmp, info) 
    464 #else 
     549#endif 
     550#if defined key_oasis3 
    465551         CALL oasis_get_freqs(id,      1, itmp, info) 
    466552#endif 
    467553         cpl_freq = itmp(1) 
     554#if defined key_oasis3mct   
     555         cpl_freq = namflddti( id )   
     556#endif  
    468557      ENDIF 
    469558      ! 
    … …  
    481570      ! 
    482571      DEALLOCATE( exfld ) 
     572      IF ( rootexchg ) DEALLOCATE ( tbuf ) 
    483573      IF (nstop == 0) THEN 
    484574         CALL oasis_terminate( nerror )          
    … …  
    489579   END SUBROUTINE cpl_finalize 
    490580 
    491 #if ! defined key_oasis3 
     581#if ! defined key_oasis3 && ! defined key_oasis3mct 
    492582 
    493583   !!---------------------------------------------------------------------- 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90

    r8058 r11286  
    5151 
    5252   SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1,   & 
    53                        px2 , py2 ) 
     53                       px2 , py2, kchoix ) 
    5454      !!---------------------------------------------------------------------- 
    5555      !!                  ***  ROUTINE repcmo  *** 
    … …  
    6767      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   px2          ! i-componante (defined at u-point) 
    6868      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   py2          ! j-componante (defined at v-point) 
    69       !!---------------------------------------------------------------------- 
    70        
    71       ! Change from geographic to stretched coordinate 
    72       ! ---------------------------------------------- 
    73       CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 ) 
    74       CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 ) 
     69      INTEGER, INTENT( IN )                       ::   kchoix       ! type of transformation    
     70                                                                    ! = 1 change from geographic to model grid.    
     71                                                                    ! =-1 change from model to geographic grid    
     72      !!----------------------------------------------------------------------   
     73         
     74      SELECT CASE (kchoix)    
     75      CASE ( 1)    
     76        ! Change from geographic to stretched coordinate    
     77        ! ----------------------------------------------    
     78         
     79        CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )    
     80        CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 )    
     81      CASE (-1)    
     82        ! Change from stretched to geographic coordinate    
     83        ! ----------------------------------------------    
     84         
     85        CALL rot_rep( pxu1, pyu1, 'U', 'ij->e',px2 )    
     86        CALL rot_rep( pxv1, pyv1, 'V', 'ij->n',py2 )    
     87      END SELECT    
    7588       
    7689   END SUBROUTINE repcmo 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r8059 r11286  
    3535   LOGICAL , PUBLIC ::   ln_blk_core    !: CORE bulk formulation 
    3636   LOGICAL , PUBLIC ::   ln_blk_mfs     !: MFS  bulk formulation 
    37 #if defined key_oasis3 
     37#if defined key_oasis3 || defined key_oasis3mct 
    3838   LOGICAL , PUBLIC ::   lk_oasis = .TRUE.  !: OASIS used 
    3939#else 
    … …  
    4242   LOGICAL , PUBLIC ::   ln_cpl         !: ocean-atmosphere coupled formulation 
    4343   LOGICAL , PUBLIC ::   ln_mixcpl      !: ocean-atmosphere forced-coupled mixed formulation 
     44   LOGICAL , PUBLIC ::   ln_wavcpl      !: ocean-wave forced-coupled mixed formulation  
     45   LOGICAL , PUBLIC ::   ll_purecpl     !: ocean-atmosphere or ocean-wave pure coupled formulation 
    4446   LOGICAL , PUBLIC ::   ln_dm2dc       !: Daily mean to Diurnal Cycle short wave (qsr) 
    4547   LOGICAL , PUBLIC ::   ln_rnf         !: runoffs / runoff mouths 
    4648   LOGICAL , PUBLIC ::   ln_ssr         !: Sea Surface restoring on SST and/or SSS       
    4749   LOGICAL , PUBLIC ::   ln_apr_dyn     !: Atmospheric pressure forcing used on dynamics (ocean & ice) 
     50   LOGICAL, PUBLIC  ::   ln_usernfmask = .false.   !  use a runoff mask file to merge data received from several models  
     51                                                   !   -> file rnfmask.nc with the float variable called rnfmask (jpi,jpj,nn_cplmodel) 
    4852   INTEGER , PUBLIC ::   nn_ice         !: flag for ice in the surface boundary condition (=0/1/2/3) 
    4953   INTEGER , PUBLIC ::   nn_isf         !: flag for isf in the surface boundary condition (=0/1/2/3/4)  
    … …  
    6468   LOGICAL , PUBLIC ::   ln_wave        !: true if some coupling with wave model 
    6569   LOGICAL , PUBLIC ::   ln_cdgw        !: true if neutral drag coefficient from wave model 
    66    LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model 
     70   LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d Stokes drift from wave model  
     71   LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used 
     72   LOGICAL , PUBLIC ::   ln_tauoc       !: true if normalized stress from wave is used   
     73   LOGICAL , PUBLIC ::   ln_tauw        !: true if ocean stress components from wave is used   
     74   LOGICAL , PUBLIC ::   ln_phioc       !: true if wave energy to ocean is used   
     75   LOGICAL , PUBLIC ::   ln_rough       !: true if wave roughness equals significant wave height  
     76   LOGICAL , PUBLIC ::   ln_zdfqiao     !: Enhanced wave vertical mixing Qiao(2010) formulation flag  
     77   INTEGER , PUBLIC ::   nn_drag        ! type of formula to calculate wind stress from wind components  
     78   INTEGER , PUBLIC ::   nn_wmix        ! type of wave breaking mixing 
    6779   ! 
    6880   LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs 
    … …  
    91103   INTEGER , PUBLIC, PARAMETER ::   jp_iam_sas  = 2      !: Multi executable configuration - SAS component 
    92104                                                         !  (internal OASIS coupling) 
     105   !!----------------------------------------------------------------------  
     106   !!           wind stress definition  
     107   !!----------------------------------------------------------------------  
     108   INTEGER, PUBLIC, PARAMETER ::   jp_ukmo  = 0        ! UKMO SHELF formulation  
     109   INTEGER, PUBLIC, PARAMETER ::   jp_std   = 1        ! standard formulation with forced or coupled drag coefficient   
     110   INTEGER, PUBLIC, PARAMETER ::   jp_const = 2        ! standard formulation with constant drag coefficient   
     111   INTEGER, PUBLIC, PARAMETER ::   jp_mcore = 3        ! momentum calculated from core forcing fields   
     112 
     113   !!----------------------------------------------------------------------  
     114   !!           Wave mixing vertical parameterization  
     115   !!----------------------------------------------------------------------  
     116   INTEGER, PUBLIC, PARAMETER ::   jp_craigbanner = 0  ! Craig and Banner formulation (original NEMO formulation -  
     117                                                       ! direct conversion of mechanical to turbulent energy)  
     118   INTEGER, PUBLIC, PARAMETER ::   jp_janssen     = 1  ! Janssen formulation - no assumption on direct energy conversion 
    93119   !!---------------------------------------------------------------------- 
    94120   !!              Ocean Surface Boundary Condition fields 
    … …  
    128154#endif 
    129155   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl) 
     156   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xrnfmask          !: coupling mask for ln_usernfmask (warning: allocated in sbcrnf) 
    130157 
    131158   !!---------------------------------------------------------------------- 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90

    r8058 r11286  
    2626    
    2727   !                                !!* namsbc_apr namelist (Atmospheric PRessure) * 
     28   LOGICAL, PUBLIC ::   cpl_mslp = .FALSE. ! Is the pressure read from coupling? 
    2829   LOGICAL, PUBLIC ::   ln_apr_obc   !: inverse barometer added to OBC ssh data  
    2930   LOGICAL, PUBLIC ::   ln_ref_apr   !: ref. pressure: global mean Patm (F) or a constant (F) 
    30    REAL(wp)        ::   rn_pref      !  reference atmospheric pressure   [N/m2] 
     31   REAL(wp),PUBLIC ::   rn_pref      !  reference atmospheric pressure   [N/m2] 
    3132 
    3233   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   ssh_ib    ! Inverse barometer now    sea surface height   [m] 
    … …  
    3435   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   apr       ! atmospheric pressure at kt                 [N/m2] 
    3536    
    36    REAL(wp) ::   tarea                ! whole domain mean masked ocean surface 
    37    REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0) 
     37   REAL(wp), PUBLIC ::   tarea                ! whole domain mean masked ocean surface 
     38   REAL(wp), PUBLIC ::   r1_grau              ! = 1.e0 / (grav * rau0) 
    3839    
    3940   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_apr   ! structure of input fields (file informations, fields read) 
    … …  
    8586         IF(lwm) WRITE ( numond, namsbc_apr ) 
    8687         ! 
    87          ALLOCATE( sf_apr(1), STAT=ierror )           !* allocate and fill sf_sst (forcing structure) with sn_sst 
    88          IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_apr: unable to allocate sf_apr structure' ) 
    89          ! 
    90          CALL fld_fill( sf_apr, (/ sn_apr /), cn_dir, 'sbc_apr', 'Atmospheric pressure ', 'namsbc_apr' ) 
    91                                 ALLOCATE( sf_apr(1)%fnow(jpi,jpj,1)   ) 
    92          IF( sn_apr%ln_tint )   ALLOCATE( sf_apr(1)%fdta(jpi,jpj,1,2) ) 
     88         IF( .NOT. cpl_mslp ) THEN  
     89            ALLOCATE( sf_apr(1), STAT=ierror )           !* allocate and fill sf_sst (forcing structure) with sn_sst  
     90            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_apr: unable to allocate sf_apr structure' )  
     91            !  
     92            CALL fld_fill( sf_apr, (/ sn_apr /), cn_dir, 'sbc_apr', 'Atmospheric pressure ', 'namsbc_apr' )  
     93                                   ALLOCATE( sf_apr(1)%fnow(jpi,jpj,1)   )  
     94            IF( sn_apr%ln_tint )   ALLOCATE( sf_apr(1)%fdta(jpi,jpj,1,2) )  
     95         ENDIF 
    9396                                ALLOCATE( ssh_ib(jpi,jpj) , ssh_ibb(jpi,jpj) ) 
    9497                                ALLOCATE( apr (jpi,jpj) ) 
    … …  
    9699         IF(lwp) THEN                                 !* control print 
    97100            WRITE(numout,*) 
    98             WRITE(numout,*) '   Namelist namsbc_apr : Atmospheric PRessure as extrenal forcing' 
     101            IF( cpl_mslp ) THEN  
     102               WRITE(numout,*) '   Namelist namsbc_apr : Atmospheric Pressure as extrenal coupling'  
     103            ELSE  
     104               WRITE(numout,*) '   Namelist namsbc_apr : Atmospheric Pressure as extrenal forcing'  
     105            ENDIF 
    99106            WRITE(numout,*) '      ref. pressure: global mean Patm (T) or a constant (F)  ln_ref_apr = ', ln_ref_apr 
    100107         ENDIF 
    … …  
    119126      ENDIF 
    120127 
    121       !                                         ! ========================== ! 
    122       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN      !    At each sbc time-step   ! 
    123          !                                      ! ===========+++============ ! 
    124          ! 
    125          IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields 
    126          ! 
    127          CALL fld_read( kt, nn_fsbc, sf_apr )               !* input Patm provided at kt + nn_fsbc/2 
    128          ! 
    129          !                                                  !* update the reference atmospheric pressure (if necessary) 
    130          IF( ln_ref_apr )   rn_pref = glob_sum( sf_apr(1)%fnow(:,:,1) * e1e2t(:,:) ) / tarea 
    131          ! 
    132          !                                                  !* Patm related forcing at kt 
    133          ssh_ib(:,:) = - ( sf_apr(1)%fnow(:,:,1) - rn_pref ) * r1_grau    ! equivalent ssh (inverse barometer) 
    134          apr   (:,:) =     sf_apr(1)%fnow(:,:,1)                        ! atmospheric pressure 
    135          ! 
    136          CALL iom_put( "ssh_ib", ssh_ib )                   !* output the inverse barometer ssh 
    137       ENDIF 
    138  
    139       !                                         ! ---------------------------------------- ! 
    140       IF( kt == nit000 ) THEN                   !   set the forcing field at nit000 - 1    ! 
    141          !                                      ! ---------------------------------------- ! 
    142          !                                            !* Restart: read in restart file 
    143          IF( ln_rstart .AND. iom_varid( numror, 'ssh_ibb', ldstop = .FALSE. ) > 0 ) THEN  
    144             IF(lwp) WRITE(numout,*) 'sbc_apr:   ssh_ibb read in the restart file' 
    145             CALL iom_get( numror, jpdom_autoglo, 'ssh_ibb', ssh_ibb )   ! before inv. barometer ssh 
     128      IF( .NOT. cpl_mslp ) THEN                    ! ========================== !  
     129         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN      !    At each sbc time-step   !  
     130            !                                      ! ===========+++============ ! 
    146131            ! 
    147          ELSE                                         !* no restart: set from nit000 values 
    148             IF(lwp) WRITE(numout,*) 'sbc_apr:   ssh_ibb set to nit000 values' 
    149             ssh_ibb(:,:) = ssh_ib(:,:) 
     132            IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields  
     133            !  
     134            CALL fld_read( kt, nn_fsbc, sf_apr )               !* input Patm provided at kt + nn_fsbc/2  
     135            !  
     136            !                                                  !* update the reference atmospheric pressure (if necessary)  
     137            IF( ln_ref_apr )   rn_pref = glob_sum( sf_apr(1)%fnow(:,:,1) * e1e2t(:,:) ) / tarea  
     138            !  
     139            !                                                  !* Patm related forcing at kt  
     140            ssh_ib(:,:) = - ( sf_apr(1)%fnow(:,:,1) - rn_pref ) * r1_grau    ! equivalent ssh (inverse barometer)  
     141            apr   (:,:) =     sf_apr(1)%fnow(:,:,1)                        ! atmospheric pressure  
     142            !  
     143            CALL iom_put( "ssh_ib", ssh_ib )                   !* output the inverse barometer ssh 
    150144         ENDIF 
    151       ENDIF 
    152       !                                         ! ---------------------------------------- ! 
    153       IF( lrst_oce ) THEN                       !      Write in the ocean restart file     ! 
    154          !                                      ! ---------------------------------------- ! 
    155          IF(lwp) WRITE(numout,*) 
    156          IF(lwp) WRITE(numout,*) 'sbc_apr : ssh_ib written in ocean restart file at it= ', kt,' date= ', ndastp 
    157          IF(lwp) WRITE(numout,*) '~~~~' 
    158          CALL iom_rstput( kt, nitrst, numrow, 'ssh_ibb' , ssh_ib ) 
     145  
     146         !                                         ! ---------------------------------------- !  
     147         IF( kt == nit000 ) THEN                   !   set the forcing field at nit000 - 1    !  
     148            !                                      ! ---------------------------------------- !  
     149            !                                            !* Restart: read in restart file  
     150            IF( ln_rstart .AND. iom_varid( numror, 'ssh_ibb', ldstop = .FALSE. ) > 0 ) THEN   
     151               IF(lwp) WRITE(numout,*) 'sbc_apr:   ssh_ibb read in the restart file'  
     152               CALL iom_get( numror, jpdom_autoglo, 'ssh_ibb', ssh_ibb )   ! before inv. barometer ssh  
     153               !  
     154            ELSE                                         !* no restart: set from nit000 values  
     155               IF(lwp) WRITE(numout,*) 'sbc_apr:   ssh_ibb set to nit000 values'  
     156               ssh_ibb(:,:) = ssh_ib(:,:)  
     157            ENDIF  
     158         ENDIF  
     159         !                                         ! ---------------------------------------- !  
     160         IF( lrst_oce ) THEN                       !      Write in the ocean restart file     !  
     161            !                                      ! ---------------------------------------- !  
     162            IF(lwp) WRITE(numout,*)  
     163            IF(lwp) WRITE(numout,*) 'sbc_apr : ssh_ib written in ocean restart file at it= ', kt,' date= ', ndastp  
     164            IF(lwp) WRITE(numout,*) '~~~~'  
     165            CALL iom_rstput( kt, nitrst, numrow, 'ssh_ibb' , ssh_ib )  
     166         ENDIF 
    159167      ENDIF 
    160168      ! 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r8058 r11286  
    319319         &               Cd, Ch, Ce, zt_zu, zq_zu ) 
    320320     
     321      IF ( ln_cdgw .AND. nn_drag==jp_std ) Cd(:,:) = cdn_wave(:,:) 
    321322      ! ... tau module, i and j component 
    322323      DO jj = 1, jpj 
    … …  
    737738 
    738739      !! Neutral coefficients at 10m: 
    739       IF( ln_cdgw ) THEN      ! wave drag case 
     740      IF( ln_cdgw .AND. nn_drag==jp_mcore ) THEN      ! wave drag case 
    740741         cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) 
    741742         ztmp0   (:,:) = cdn_wave(:,:) 
    … …  
    783784         END IF 
    784785        
    785          IF( ln_cdgw ) THEN      ! surface wave case 
     786         IF( ln_cdgw .AND. nn_drag==jp_mcore ) THEN      ! surface wave case 
    786787            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u )  
    787788            Cd      = sqrt_Cd * sqrt_Cd 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90

    r8058 r11286  
    2424   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2525   USE prtctl          ! Print control 
    26    USE sbcwave,ONLY : cdn_wave !wave module 
    2726 
    2827   IMPLICIT NONE 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r8058 r11286  
    2323   USE sbcapr 
    2424   USE sbcdcy          ! surface boundary condition: diurnal cycle 
     25   USE sbcwave         ! surface boundary condition: waves 
    2526   USE phycst          ! physical constants 
    2627#if defined key_lim3 
    … …  
    4142   USE timing          ! Timing 
    4243   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     44#if defined key_oasis3 || defined key_oasis3mct    
     45   USE mod_oasis                    ! OASIS3-MCT module    
     46#endif   
    4347   USE eosbn2 
    4448   USE sbcrnf   , ONLY : l_rnfcpl 
    … …  
    105109   INTEGER, PARAMETER ::   jpr_e3t1st = 41            ! first T level thickness  
    106110   INTEGER, PARAMETER ::   jpr_fraqsr = 42            ! fraction of solar net radiation absorbed in the first ocean level 
    107    INTEGER, PARAMETER ::   jprcv      = 42            ! total number of fields received 
     111   INTEGER, PARAMETER ::   jpr_mslp   = 43            ! mean sea level pressure    
     112   INTEGER, PARAMETER ::   jpr_hsig   = 44            ! Hsig    
     113   INTEGER, PARAMETER ::   jpr_phioc  = 45            ! Wave=>ocean energy flux    
     114   INTEGER, PARAMETER ::   jpr_sdrftx = 46            ! Stokes drift on grid 1    
     115   INTEGER, PARAMETER ::   jpr_sdrfty = 47            ! Stokes drift on grid 2    
     116   INTEGER, PARAMETER ::   jpr_wper   = 48            ! Mean wave period   
     117   INTEGER, PARAMETER ::   jpr_wnum   = 49            ! Mean wavenumber   
     118   INTEGER, PARAMETER ::   jpr_tauoc  = 50            ! Stress fraction adsorbed by waves 
     119   INTEGER, PARAMETER ::   jpr_wdrag  = 51            ! Neutral surface drag coefficient   
     120   INTEGER, PARAMETER ::   jpr_wfreq  = 52            ! Wave peak frequency   
     121   INTEGER, PARAMETER ::   jpr_tauwx  = 53            ! x component of the ocean stress from waves  
     122   INTEGER, PARAMETER ::   jpr_tauwy  = 54            ! y component of the ocean stress from waves  
     123   INTEGER, PARAMETER ::   jprcv      = 54            ! total number of fields received 
    108124 
    109125   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere 
    … …  
    135151   INTEGER, PARAMETER ::   jps_e3t1st = 27            ! first level depth (vvl) 
    136152   INTEGER, PARAMETER ::   jps_fraqsr = 28            ! fraction of solar net radiation absorbed in the first ocean level 
    137    INTEGER, PARAMETER ::   jpsnd      = 28            ! total number of fields sended 
     153   INTEGER, PARAMETER ::   jps_ficet  = 29            ! total ice fraction     
     154   INTEGER, PARAMETER ::   jps_ocxw   = 30            ! currents on grid 1     
     155   INTEGER, PARAMETER ::   jps_ocyw   = 31            ! currents on grid 2   
     156   INTEGER, PARAMETER ::   jps_wlev   = 32            ! water level    
     157   INTEGER, PARAMETER ::   jpsnd      = 32            ! total number of fields sent 
    138158 
    139159   !                                                         !!** namelist namsbc_cpl ** 
    … …  
    149169   ! Received from the atmosphere                     ! 
    150170   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf 
    151    TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                         
     171   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp                              
     172   ! Send to waves    
     173   TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev    
     174   ! Received from waves    
     175   TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrft,sn_rcv_wper, &  
     176                    sn_rcv_wfreq,sn_rcv_wnum,sn_rcv_tauoc,sn_rcv_tauw, &  
     177                    sn_rcv_wdrag 
    152178   ! Other namelist parameters                        ! 
    153179   INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 
    … …  
    161187 
    162188   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky) 
    163  
     189     
    164190   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument 
    165191 
    … …  
    179205      !!             ***  FUNCTION sbc_cpl_alloc  *** 
    180206      !!---------------------------------------------------------------------- 
    181       INTEGER :: ierr(3) 
     207      INTEGER :: ierr(4) 
    182208      !!---------------------------------------------------------------------- 
    183209      ierr(:) = 0 
    … …  
    188214      ALLOCATE( a_i(jpi,jpj,1) , STAT=ierr(2) )  ! used in sbcice_if.F90 (done here as there is no sbc_ice_if_init) 
    189215#endif 
    190       ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) ) 
    191       ! 
     216!      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) ) 
     217      ! Hardwire three models as nn_cplmodel has not been read in from the namelist yet.    
     218      ALLOCATE( xcplmask(jpi,jpj,0:3) , STAT=ierr(3) ) 
     219      ! 
     220      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) ) 
     221 
    192222      sbc_cpl_alloc = MAXVAL( ierr ) 
    193223      IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc ) 
    … …  
    216246      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos 
    217247      !! 
    218       NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,      & 
    219          &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,      & 
    220          &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx,   & 
    221          &                  sn_rcv_co2 , nn_cplmodel  , ln_usecplmask 
     248      NAMELIST/namsbc_cpl/  sn_snd_temp , sn_snd_alb  , sn_snd_thick , sn_snd_crt   , sn_snd_co2,      &    
     249         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,      &    
     250         &                  sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev  , sn_rcv_hsig  , sn_rcv_phioc ,   &    
     251         &                  sn_rcv_sdrft, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_wfreq, sn_rcv_tauoc,    &  
     252         &                  sn_rcv_wdrag, sn_rcv_qns   , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal ,     &  
     253         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp , sn_rcv_tauw ,                  &  
     254         &                  nn_cplmodel, ln_usecplmask, ln_usernfmask 
    222255      !!--------------------------------------------------------------------- 
    223256      ! 
    … …  
    260293         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 
    261294         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')' 
     295         WRITE(numout,*)'      mean sea level pressure         = ', TRIM(sn_rcv_mslp%cldes  ), ' (', TRIM(sn_rcv_mslp%clcat  ), ')' 
     296         WRITE(numout,*)'      significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')'    
     297         WRITE(numout,*)'      wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')'    
     298         WRITE(numout,*)'      Surface Stokes drift u,v        = ', TRIM(sn_rcv_sdrft%cldes ), ' (', TRIM(sn_rcv_sdrft%clcat ), ')' 
     299         WRITE(numout,*)'      Mean wave period                = ', TRIM(sn_rcv_wper%cldes  ), ' (', TRIM(sn_rcv_wper%clcat  ), ')'    
     300         WRITE(numout,*)'      Mean wave number                = ', TRIM(sn_rcv_wnum%cldes  ), ' (', TRIM(sn_rcv_wnum%clcat  ), ')'    
     301         WRITE(numout,*)'      Wave peak frequency             = ', TRIM(sn_rcv_wfreq%cldes ), ' (', TRIM(sn_rcv_wfreq%clcat ), ')'    
     302         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_tauoc%cldes ), ' (', TRIM(sn_rcv_tauoc%clcat ), ')'    
     303         WRITE(numout,*)'      Stress components by waves      = ', TRIM(sn_rcv_tauw%cldes  ), ' (', TRIM(sn_rcv_tauw%clcat  ), ')' 
     304         WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 
    262305         WRITE(numout,*)'  sent fields (multiple ice categories)' 
    263306         WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')' 
    264307         WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')' 
    265308         WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')' 
     309         WRITE(numout,*)'      total ice fraction              = ', TRIM(sn_snd_ifrac%cldes ), ' (', TRIM(sn_snd_ifrac%clcat ), ')' 
    266310         WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')' 
    267311         WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref  
    … …  
    269313         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd 
    270314         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')' 
     315         WRITE(numout,*)'      water level                     = ', TRIM(sn_snd_wlev%cldes  ), ' (', TRIM(sn_snd_wlev%clcat  ), ')'    
     316         WRITE(numout,*)'      surface current to waves        = ', TRIM(sn_snd_crtw%cldes  ), ' (', TRIM(sn_snd_crtw%clcat  ), ')'    
     317         WRITE(numout,*)'                      - referential   = ', sn_snd_crtw%clvref    
     318         WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor    
     319         WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd 
    271320         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel 
    272321         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask 
     322         WRITE(numout,*)'  ln_usernfmask                       = ', ln_usernfmask 
    273323      ENDIF 
    274324 
    275325      !                                   ! allocate sbccpl arrays 
    276       IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' ) 
     326      !IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' ) 
    277327      
    278328      ! ================================ ! 
    … …  
    307357      !  
    308358      ! Vectors: change of sign at north fold ONLY if on the local grid 
     359      IF( TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM(sn_rcv_tau%cldes ) == 'oce and ice') THEN ! avoid working with the atmospheric fields if they are not coupled 
    309360      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1. 
    310361       
    … …  
    337388         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point 
    338389         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point 
    339          srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2 
     390         !srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2   
     391         ! Currently needed for HadGEM3 - but shouldn't affect anyone else for the moment   
     392         srcv(jpr_otx1)%laction = .TRUE.    
     393         srcv(jpr_oty1)%laction = .TRUE.   
     394         !   
    340395         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only 
    341396      CASE( 'T,I' )  
    … …  
    374429         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp. 
    375430      ENDIF 
     431      ENDIF 
    376432        
    377433      !                                                      ! ------------------------- ! 
    … …  
    403459      IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN 
    404460         srcv(jpr_rnf)%laction = .TRUE. 
    405          l_rnfcpl              = .TRUE.                      ! -> no need to read runoffs in sbcrnf 
     461         l_rnfcpl              = .NOT. ln_usernfmask         ! -> no need to read runoffs in sbcrnf 
    406462         ln_rnf                = nn_components /= jp_iam_sas ! -> force to go through sbcrnf if not sas 
    407463         IF(lwp) WRITE(numout,*) 
    … …  
    470526      !                                                      ! ------------------------- ! 
    471527      srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE. 
     528 
     529      !                                                      ! ------------------------- !    
     530      !                                                      ! Mean Sea Level Pressure   !    
     531      !                                                      ! ------------------------- !    
     532      srcv(jpr_mslp)%clname = 'O_MSLP'  
     533      IF( TRIM(sn_rcv_mslp%cldes  ) == 'coupled' ) THEN  
     534         srcv(jpr_mslp)%laction = .TRUE.  
     535         cpl_mslp = .TRUE.  
     536      ENDIF 
     537 
    472538      !                                                      ! ------------------------- ! 
    473539      !                                                      !   topmelt and botmelt     !    
    … …  
    483549         srcv(jpr_topm:jpr_botm)%laction = .TRUE. 
    484550      ENDIF 
     551      !                                                      ! ------------------------- !   
     552      !                                                      !      Wave breaking        !       
     553      !                                                      ! ------------------------- !    
     554      srcv(jpr_hsig)%clname  = 'O_Hsigwa'    ! significant wave height   
     555      IF( TRIM(sn_rcv_hsig%cldes  ) == 'coupled' )  THEN   
     556         srcv(jpr_hsig)%laction = .TRUE.   
     557         cpl_hsig = .TRUE.   
     558      ENDIF   
     559      srcv(jpr_phioc)%clname = 'O_PhiOce'    ! wave to ocean energy   
     560      IF( TRIM(sn_rcv_phioc%cldes ) == 'coupled' )  THEN   
     561         srcv(jpr_phioc)%laction = .TRUE.   
     562         cpl_phioc = .TRUE.   
     563      ENDIF   
     564      srcv(jpr_sdrftx)%clname = 'O_Sdrfx'    ! Stokes drift in the u direction   
     565      srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction   
     566      IF( TRIM(sn_rcv_sdrft%cldes ) == 'coupled' )  THEN 
     567         srcv(jpr_sdrftx)%laction = .TRUE.   
     568         srcv(jpr_sdrfty)%laction = .TRUE.   
     569         cpl_sdrft = .TRUE. 
     570      ENDIF   
     571      srcv(jpr_wper)%clname = 'O_WPer'       ! mean wave period   
     572      IF( TRIM(sn_rcv_wper%cldes  ) == 'coupled' )  THEN   
     573         srcv(jpr_wper)%laction = .TRUE.   
     574         cpl_wper = .TRUE.   
     575      ENDIF   
     576      srcv(jpr_wfreq)%clname = 'O_WFreq'     ! wave peak frequency   
     577      IF( TRIM(sn_rcv_wfreq%cldes ) == 'coupled' )  THEN   
     578         srcv(jpr_wfreq)%laction = .TRUE.   
     579         cpl_wfreq = .TRUE.   
     580      ENDIF 
     581      srcv(jpr_wnum)%clname = 'O_WNum'       ! mean wave number   
     582      IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' )  THEN   
     583         srcv(jpr_wnum)%laction = .TRUE.   
     584         cpl_wnum = .TRUE.   
     585      ENDIF   
     586      srcv(jpr_tauoc)%clname = 'O_TauOce'     ! stress fraction adsorbed by the wave   
     587      IF( TRIM(sn_rcv_tauoc%cldes ) == 'coupled' )  THEN   
     588         srcv(jpr_tauoc)%laction = .TRUE.   
     589         cpl_tauoc = .TRUE.   
     590      ENDIF   
     591      srcv(jpr_tauwx)%clname = 'O_Tauwx'      ! ocean stress from wave in the x direction  
     592      srcv(jpr_tauwy)%clname = 'O_Tauwy'      ! ocean stress from wave in the y direction  
     593      IF( TRIM(sn_rcv_tauw%cldes ) == 'coupled' )  THEN   
     594         srcv(jpr_tauwx)%laction = .TRUE.   
     595         srcv(jpr_tauwy)%laction = .TRUE.   
     596         cpl_tauw = .TRUE. 
     597      ENDIF   
     598      srcv(jpr_wdrag)%clname = 'O_WDrag'     ! neutral surface drag coefficient   
     599      IF( TRIM(sn_rcv_wdrag%cldes ) == 'coupled' )  THEN   
     600         srcv(jpr_wdrag)%laction = .TRUE.   
     601         cpl_wdrag = .TRUE.   
     602      ENDIF   
     603      !  
     604      IF( srcv(jpr_tauoc)%laction .AND. srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction ) &  
     605            CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', &  
     606                                     '(sn_rcv_tauoc=coupled and sn_rcv_tauw=coupled)' )  
     607      ! 
     608      !    
    485609      !                                                      ! ------------------------------- ! 
    486610      !                                                      !   OPA-SAS coupling - rcv by opa !    
    … …  
    637761      !                                                      ! ------------------------- ! 
    638762      ssnd(jps_fice)%clname = 'OIceFrc' 
     763      ssnd(jps_ficet)%clname = 'OIceFrcT' 
    639764      ssnd(jps_hice)%clname = 'OIceTck' 
    640765      ssnd(jps_hsnw)%clname = 'OSnwTck' 
    … …  
    645770      ENDIF 
    646771       
     772      IF (TRIM( sn_snd_ifrac%cldes )  == 'coupled') ssnd(jps_ficet)%laction = .TRUE. 
     773 
    647774      SELECT CASE ( TRIM( sn_snd_thick%cldes ) ) 
    648775      CASE( 'none'         )       ! nothing to do 
    … …  
    665792      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1' 
    666793      ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1' 
     794      ssnd(jps_ocxw)%clname = 'O_OCurxw'    
     795      ssnd(jps_ocyw)%clname = 'O_OCuryw' 
    667796      ! 
    668797      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold 
    … …  
    685814      END SELECT 
    686815 
     816      ssnd(jps_ocxw:jps_ocyw)%nsgn = -1.   ! vectors: change of the sign at the north fold    
     817               
     818      IF( sn_snd_crtw%clvgrd == 'U,V' ) THEN    
     819         ssnd(jps_ocxw)%clgrid = 'U' ; ssnd(jps_ocyw)%clgrid = 'V'    
     820      ELSE IF( sn_snd_crtw%clvgrd /= 'T' ) THEN    
     821         CALL ctl_stop( 'sn_snd_crtw%clvgrd must be equal to T' )    
     822      ENDIF    
     823      IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) ssnd(jps_ocxw:jps_ocyw)%nsgn = 1.    
     824      SELECT CASE( TRIM( sn_snd_crtw%cldes ) )    
     825         CASE( 'none'                 )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .FALSE.    
     826         CASE( 'oce only'             )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .TRUE.    
     827         CASE( 'weighted oce and ice' )   !   nothing to do    
     828         CASE( 'mixed oce-ice'        )   ; ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.    
     829         CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crtw%cldes' )    
     830      END SELECT 
     831 
    687832      !                                                      ! ------------------------- ! 
    688833      !                                                      !          CO2 flux         ! 
    689834      !                                                      ! ------------------------- ! 
    690835      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE. 
     836 
     837      !                                                      ! ------------------------- !    
     838      !                                                      !     Sea surface height    !    
     839      !                                                      ! ------------------------- !    
     840      ssnd(jps_wlev)%clname = 'O_Wlevel' ;  IF( TRIM(sn_snd_wlev%cldes) == 'coupled' )   ssnd(jps_wlev)%laction = .TRUE. 
    691841 
    692842      !                                                      ! ------------------------------- ! 
    … …  
    783933      IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   & 
    784934         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' ) 
    785       ncpl_qsr_freq = 86400 / ncpl_qsr_freq 
     935      IF( ln_dm2dc .AND. ln_cpl ) ncpl_qsr_freq = 86400 / ncpl_qsr_freq 
    786936 
    787937      CALL wrk_dealloc( jpi,jpj, zacs, zaos ) 
    … …  
    837987      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case) 
    838988      !!---------------------------------------------------------------------- 
     989      USE sbcflx ,  ONLY : ln_shelf_flx  
     990      USE sbcssm ,  ONLY : sbc_ssm_cpl  
     991      USE lib_fortran     ! distributed memory computing library 
     992 
    839993      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index 
    840994      INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation  
    … …  
    845999      INTEGER  ::   ji, jj, jn             ! dummy loop indices 
    8461000      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000) 
     1001      INTEGER  ::   ikchoix 
    8471002      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars      
    8481003      REAL(wp) ::   zcoef                  ! temporary scalar 
    … …  
    8501005      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient 
    8511006      REAL(wp) ::   zzx, zzy               ! temporary variables 
    852       REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty, zmsk, zemp, zqns, zqsr 
     1007      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty, ztx2, zty2, zmsk, zemp, zqns, zqsr 
    8531008      !!---------------------------------------------------------------------- 
    8541009      ! 
    8551010      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv') 
    8561011      ! 
    857       CALL wrk_alloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr ) 
     1012      CALL wrk_alloc( jpi,jpj, ztx, zty, ztx2, zty2, zmsk, zemp, zqns, zqsr ) 
    8581013      ! 
    8591014      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0) 
    … …  
    8931048            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid 
    8941049               !                                                       ! (geographical to local grid -> rotate the components) 
    895                CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )    
    896                IF( srcv(jpr_otx2)%laction ) THEN 
    897                   CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )    
    898                ELSE   
    899                   CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )   
     1050               IF( srcv(jpr_otx1)%clgrid == 'U' .AND. (.NOT. srcv(jpr_otx2)%laction) ) THEN    
     1051                  ! Temporary code for HadGEM3 - will be removed eventually.    
     1052                  ! Only applies when we have only taux on U grid and tauy on V grid    
     1053                  DO jj=2,jpjm1    
     1054                     DO ji=2,jpim1    
     1055                        ztx(ji,jj)=0.25*vmask(ji,jj,1) &    
     1056                           *(frcv(jpr_otx1)%z3(ji,jj,1)+frcv(jpr_otx1)%z3(ji-1,jj,1)    &    
     1057                           +frcv(jpr_otx1)%z3(ji,jj+1,1)+frcv(jpr_otx1)%z3(ji-1,jj+1,1))    
     1058                        zty(ji,jj)=0.25*umask(ji,jj,1) &    
     1059                           *(frcv(jpr_oty1)%z3(ji,jj,1)+frcv(jpr_oty1)%z3(ji+1,jj,1)    &    
     1060                           +frcv(jpr_oty1)%z3(ji,jj-1,1)+frcv(jpr_oty1)%z3(ji+1,jj-1,1))    
     1061                     ENDDO    
     1062                  ENDDO    
     1063                               
     1064                  ikchoix = 1    
     1065                  CALL repcmo(frcv(jpr_otx1)%z3(:,:,1),zty,ztx,frcv(jpr_oty1)%z3(:,:,1),ztx2,zty2,ikchoix)    
     1066                  CALL lbc_lnk (ztx2,'U', -1. )    
     1067                  CALL lbc_lnk (zty2,'V', -1. )    
     1068                  frcv(jpr_otx1)%z3(:,:,1)=ztx2(:,:)    
     1069                  frcv(jpr_oty1)%z3(:,:,1)=zty2(:,:)    
     1070               ELSE    
     1071                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )       
     1072                  frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid    
     1073                  IF( srcv(jpr_otx2)%laction ) THEN    
     1074                     CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )       
     1075                  ELSE    
     1076                     CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )     
     1077                  ENDIF    
     1078                  frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid 
    9001079               ENDIF 
    901                frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid 
    902                frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid 
    9031080            ENDIF 
    9041081            !                               
    … …  
    9191096      ELSE                                                   !   No dynamical coupling   ! 
    9201097         !                                                   ! ========================= ! 
     1098         ! it is possible that the momentum is calculated from the winds (ln_shelf_flx) and a coupled drag coefficient  
     1099         IF( srcv(jpr_wdrag)%laction .AND. ln_shelf_flx .AND. ln_cdgw .AND. nn_drag == jp_std ) THEN  
     1100            DO jj = 1, jpjm1  
     1101               DO ji = 1, jpim1  
     1102                  ! here utau and vtau should contain the wind components as read from the forcing files,  
     1103                  ! and the wind module should already be properly calculated  
     1104                  frcv(jpr_otx1)%z3(ji,jj,1) = zrhoa * 0.5 * ( frcv(jpr_wdrag)%z3(ji,jj,1) + frcv(jpr_wdrag)%z3(ji+1,jj,1) ) * &  
     1105                                                                         utau(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) )  
     1106                  frcv(jpr_oty1)%z3(ji,jj,1) = zrhoa * 0.5 * ( frcv(jpr_wdrag)%z3(ji,jj,1) + frcv(jpr_wdrag)%z3(ji,jj+1,1) ) * &  
     1107                                                                         vtau(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) )  
     1108                  utau(ji,jj) = frcv(jpr_otx1)%z3(ji,jj,1)  
     1109                  vtau(ji,jj) = frcv(jpr_oty1)%z3(ji,jj,1)  
     1110               END DO  
     1111            END DO  
     1112            CALL lbc_lnk_multi( frcv(jpr_otx1)%z3(:,:,1), 'U', -1. , frcv(jpr_oty1)%z3(:,:,1), 'V', -1. , &  
     1113                                                             utau(:,:), 'U', -1. , vtau(:,:), 'V',  -1. )  
     1114            llnewtx = .TRUE.  
     1115         ELSE 
    9211116         frcv(jpr_otx1)%z3(:,:,1) = 0.e0                               ! here simply set to zero  
    9221117         frcv(jpr_oty1)%z3(:,:,1) = 0.e0                               ! an external read in a file can be added instead 
    9231118         llnewtx = .TRUE. 
     1119         ENDIF 
    9241120         ! 
    9251121      ENDIF 
    … …  
    9411137            END DO 
    9421138            CALL lbc_lnk( frcv(jpr_taum)%z3(:,:,1), 'T', 1. ) 
     1139            IF( .NOT. srcv(jpr_otx1)%laction .AND. srcv(jpr_wdrag)%laction .AND. &  
     1140                                ln_shelf_flx .AND. ln_cdgw .AND. nn_drag == jp_std ) &  
     1141               taum(:,:) = frcv(jpr_taum)%z3(:,:,1) 
    9431142            llnewtau = .TRUE. 
    9441143         ELSE 
    … …  
    9551154      !                                                      ! ========================= ! 
    9561155      !                                                      !      10 m wind speed      !   (wndm) 
     1156      !                                                      !   include wave drag coef  !   (wndm) 
    9571157      !                                                      ! ========================= ! 
    9581158      ! 
    … …  
    9651165!CDIR NOVERRCHK 
    9661166               DO ji = 1, jpi  
     1167                  IF( ln_shelf_flx ) THEN   ! the 10 wind module is properly calculated before if ln_shelf_flx  
     1168                     frcv(jpr_w10m)%z3(ji,jj,1) = wndm(ji,jj)  
     1169                  ELSE 
    9671170                  frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef ) 
     1171                  ENDIF 
    9681172               END DO 
    9691173            END DO 
    … …  
    9751179      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN 
    9761180         ! 
     1181         ! if ln_wavcpl, the fields already contain the right information from forcing even if not ln_mixcpl 
    9771182         IF( ln_mixcpl ) THEN 
    978             utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:) 
    979             vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:) 
    980             taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:) 
    981             wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:) 
    982          ELSE 
     1183            IF( srcv(jpr_otx1)%laction ) THEN  
     1184               utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:)  
     1185               vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:)  
     1186            ENDIF  
     1187            IF( srcv(jpr_taum)%laction )   &  
     1188               taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:)  
     1189            IF( srcv(jpr_w10m)%laction )   &  
     1190               wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:)  
     1191         ELSE IF( ll_purecpl ) THEN 
    9831192            utau(:,:) = frcv(jpr_otx1)%z3(:,:,1) 
    9841193            vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1) 
    … …  
    9961205      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1) 
    9971206#endif 
     1207      !    
     1208      !                                                      ! ========================= !    
     1209      !                                                      ! Mean Sea Level Pressure   !   (taum)    
     1210      !                                                      ! ========================= !    
     1211      !    
     1212      IF( srcv(jpr_mslp)%laction ) THEN                    ! UKMO SHELF effect of atmospheric pressure on SSH    
     1213          IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields    
     1214        
     1215          !                                                  !* update the reference atmospheric pressure (if necessary)  
     1216          IF( ln_ref_apr )  rn_pref = glob_sum( frcv(jpr_mslp)%z3(:,:,1) * e1e2t(:,:) ) / tarea  
     1217         
     1218          ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rn_pref ) * r1_grau    ! equivalent ssh (inverse barometer)    
     1219          apr   (:,:) =     frcv(jpr_mslp)%z3(:,:,1) !atmospheric pressure    
     1220          !  
     1221          CALL iom_put( "ssh_ib", ssh_ib )                                    !* output the inverse barometer ssh 
     1222        
     1223          !                                         ! ---------------------------------------- !  
     1224          IF( kt == nit000 ) THEN                   !   set the forcing field at nit000 - 1    !  
     1225             !                                      ! ---------------------------------------- !  
     1226             !* Restart: read in restart file  
     1227             IF( ln_rstart .AND. iom_varid( numror, 'ssh_ibb', ldstop = .FALSE. ) > 0 ) THEN  
     1228                IF(lwp) WRITE(numout,*) 'sbc_cpl:   ssh_ibb read in the restart file'  
     1229                CALL iom_get( numror, jpdom_autoglo, 'ssh_ibb', ssh_ibb )   ! before inv. barometer ssh  
     1230             ELSE                                         !* no restart: set from nit000 values  
     1231                IF(lwp) WRITE(numout,*) 'sbc_cpl:   ssh_ibb set to nit000 values'  
     1232                ssh_ibb(:,:) = ssh_ib(:,:)  
     1233             ENDIF  
     1234          ENDIF  
     1235          !                                         ! ---------------------------------------- !  
     1236          IF( lrst_oce ) THEN                       !      Write in the ocean restart file     !  
     1237             !                                      ! ---------------------------------------- !  
     1238             IF(lwp) WRITE(numout,*)  
     1239             IF(lwp) WRITE(numout,*) 'sbc_cpl : ssh_ib written in ocean restart file at it= ', kt,' date= ', ndastp  
     1240             IF(lwp) WRITE(numout,*) '~~~~'  
     1241             CALL iom_rstput( kt, nitrst, numrow, 'ssh_ibb' , ssh_ib )  
     1242          ENDIF  
     1243         
     1244          ! Update mean ssh  
     1245          CALL sbc_ssm_cpl( kt ) 
     1246      END IF    
     1247      !   
     1248      IF( ln_sdw ) THEN  ! Stokes Drift correction activated   
     1249      !                                                      ! ========================= !    
     1250      !                                                      !       Stokes drift u,v    !   
     1251      !                                                      ! ========================= !    
     1252      IF( srcv(jpr_sdrftx)%laction .AND. srcv(jpr_sdrfty)%laction ) THEN  
     1253                                     ut0sd(:,:) = frcv(jpr_sdrftx)%z3(:,:,1)   
     1254                                     vt0sd(:,:) = frcv(jpr_sdrfty)%z3(:,:,1)   
     1255      ENDIF 
     1256      !   
     1257      !                                                      ! ========================= !    
     1258      !                                                      !      Wave mean period     !   
     1259      !                                                      ! ========================= !    
     1260         IF( srcv(jpr_wper)%laction ) wmp(:,:) = frcv(jpr_wper)%z3(:,:,1)   
     1261      !   
     1262      !                                                      ! ========================= !    
     1263      !                                                      !  Significant wave height  !   
     1264      !                                                      ! ========================= !    
     1265         IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1)   
     1266      !   
     1267      !                                                      ! ========================= !    
     1268      !                                                      !    Wave peak frequency    !   
     1269      !                                                      ! ========================= !    
     1270         IF( srcv(jpr_wfreq)%laction ) wfreq(:,:) = frcv(jpr_wfreq)%z3(:,:,1)   
     1271      !  
     1272      !                                                      ! ========================= ! 
     1273      !                                                      !    Vertical mixing Qiao   !   
     1274      !                                                      ! ========================= !    
     1275         IF( srcv(jpr_wnum)%laction .AND. ln_zdfqiao ) wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1)   
     1276        
     1277         ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode   
     1278         IF( (srcv(jpr_sdrftx)%laction .AND. srcv(jpr_sdrfty)%laction) .OR. srcv(jpr_wper)%laction &   
     1279                                        .OR. srcv(jpr_hsig)%laction   .OR. srcv(jpr_wfreq)%laction) & 
     1280            CALL sbc_stokes()   
     1281      ENDIF   
     1282      !                                                      ! ========================= !    
     1283      !                                                      ! Stress adsorbed by waves  !   
     1284      !                                                      ! ========================= !    
     1285      IF( srcv(jpr_tauoc)%laction .AND. ln_tauoc ) THEN  
     1286         tauoc_wave(:,:) = frcv(jpr_tauoc)%z3(:,:,1)  
     1287         ! cap the value of tauoc  
     1288         WHERE(tauoc_wave <   0.0 ) tauoc_wave = 1.0  
     1289         WHERE(tauoc_wave > 100.0 ) tauoc_wave = 1.0  
     1290      ENDIF  
     1291      !                                                      ! ========================= !    
     1292      !                                                      ! Stress component by waves !   
     1293      !                                                      ! ========================= !    
     1294      IF( srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction .AND. ln_tauw ) THEN  
     1295         tauw_x(:,:) = frcv(jpr_tauwx)%z3(:,:,1)  
     1296         tauw_y(:,:) = frcv(jpr_tauwy)%z3(:,:,1)  
     1297         ! cap the value of tauoc  
     1298         WHERE(tauw_x < -100.0 ) tauw_x = 0.0  
     1299         WHERE(tauw_x >  100.0 ) tauw_x = 0.0  
     1300         WHERE(tauw_y < -100.0 ) tauw_y = 0.0  
     1301         WHERE(tauw_y >  100.0 ) tauw_y = 0.0  
     1302      ENDIF 
     1303        
     1304      !                                                      ! ========================= !    
     1305      !                                                      !   Wave to ocean energy    ! 
     1306      !                                                      ! ========================= !    
     1307      IF( srcv(jpr_phioc)%laction .AND. ln_phioc ) THEN  
     1308         rn_crban(:,:) = 29.0 * frcv(jpr_phioc)%z3(:,:,1)  
     1309         WHERE( rn_crban <    0.0 ) rn_crban = 0.0  
     1310         WHERE( rn_crban > 1000.0 ) rn_crban = 1000.0  
     1311      ENDIF  
    9981312 
    9991313      !  Fields received by SAS when OASIS coupling 
    … …  
    10671381               CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' ) 
    10681382            END SELECT 
    1069          ELSE 
     1383         ELSE IF( ll_purecpl ) THEN 
    10701384            zemp(:,:) = 0._wp 
    10711385         ENDIF 
    … …  
    10751389         IF( srcv(jpr_cal)%laction )     zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1) 
    10761390          
    1077          IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:) 
    1078          ELSE                   ;   emp(:,:) =                              zemp(:,:) 
     1391         IF( ln_mixcpl .AND. ( srcv(jpr_oemp)%laction .OR. srcv(jpr_rain)%laction )) THEN  
     1392                                         emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:)  
     1393         ELSE IF( ll_purecpl ) THEN  ;   emp(:,:) = zemp(:,:) 
    10791394         ENDIF 
    10801395         ! 
    … …  
    10911406            ENDIF 
    10921407         ENDIF 
    1093          IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:) 
    1094          ELSE                   ;   qns(:,:) =                              zqns(:,:) 
     1408         IF( ln_mixcpl .AND. ( srcv(jpr_qnsoce)%laction .OR. srcv(jpr_qnsmix)%laction )) THEN  
     1409                                          qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)  
     1410         ELSE IF( ll_purecpl ) THEN   ;   qns(:,:) = zqns(:,:) 
    10951411         ENDIF 
    10961412 
    … …  
    11011417         ENDIF 
    11021418         IF( ln_dm2dc .AND. ln_cpl )   zqsr(:,:) = sbc_dcy( zqsr )   ! modify qsr to include the diurnal cycle 
    1103          IF( ln_mixcpl ) THEN   ;   qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:) 
    1104          ELSE                   ;   qsr(:,:) =                              zqsr(:,:) 
     1419         IF( ln_mixcpl .AND. ( srcv(jpr_qsroce)%laction .OR. srcv(jpr_qsrmix)%laction )) THEN  
     1420                                          qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:)  
     1421         ELSE IF( ll_purecpl ) THEN   ;   qsr(:,:) = zqsr(:,:) 
    11051422         ENDIF 
    11061423         ! 
    … …  
    11131430      ENDIF 
    11141431      ! 
    1115       CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr ) 
     1432      CALL wrk_dealloc( jpi,jpj, ztx, zty, ztx2, zty2, zmsk, zemp, zqns, zqsr ) 
    11161433      ! 
    11171434      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv') 
    … …  
    17082025      ! 
    17092026      INTEGER ::   ji, jj, jl   ! dummy loop indices 
     2027      INTEGER ::   ikchoix 
    17102028      INTEGER ::   isec, info   ! local integer 
    17112029      REAL(wp) ::   zumax, zvmax 
    … …  
    17362054            ! 
    17372055            SELECT CASE( sn_snd_temp%cldes) 
    1738             CASE( 'oce only'             )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0 
     2056            CASE( 'oce only'             )   ;   ztmp1(:,:) = (ztmp1(:,:) + rt0) * tmask(:,:,1) 
    17392057            CASE( 'oce and ice'          )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0 
    17402058               SELECT CASE( sn_snd_temp%clcat ) 
    … …  
    17652083                  ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl) 
    17662084               ENDDO 
     2085            CASE( 'none'         )       ! nothing to do 
    17672086            CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' ) 
    17682087            END SELECT 
    … …  
    18892208         !                                                  j+1   j     -----V---F 
    18902209         ! surface velocity always sent from T point                     !       | 
    1891          !                                                        j      |   T   U 
     2210         ! [except for HadGEM3]                                   j      |   T   U 
    18922211         !                                                               |       | 
    18932212         !                                                   j    j-1   -I-------| 
    … …  
    19012220            SELECT CASE( TRIM( sn_snd_crt%cldes ) ) 
    19022221            CASE( 'oce only'             )      ! C-grid ==> T 
    1903                DO jj = 2, jpjm1 
    1904                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    1905                      zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) ) 
    1906                      zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) )  
    1907                   END DO 
    1908                END DO 
     2222               IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN    
     2223                  DO jj = 2, jpjm1    
     2224                     DO ji = fs_2, fs_jpim1   ! vector opt.    
     2225                        zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )    
     2226                        zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji,jj-1,1) )     
     2227                     END DO    
     2228                  END DO    
     2229               ELSE    
     2230                  ! Temporarily Changed for UKV    
     2231                  DO jj = 2, jpjm1    
     2232                     DO ji = 2, jpim1    
     2233                        zotx1(ji,jj) = un(ji,jj,1)    
     2234                        zoty1(ji,jj) = vn(ji,jj,1)    
     2235                     END DO    
     2236                  END DO    
     2237               ENDIF  
    19092238            CASE( 'weighted oce and ice' )    
    19102239               SELECT CASE ( cp_ice_msh ) 
    … …  
    19302259                  END DO 
    19312260               CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T 
    1932                   DO jj = 2, jpjm1 
    1933                      DO ji = 2, jpim1   ! NO vector opt. 
    1934                         zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   
    1935                         zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   
    1936                         zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     & 
    1937                            &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
    1938                         zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     & 
    1939                            &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
    1940                      END DO 
    1941                   END DO 
     2261                  IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN    
     2262                     DO jj = 2, jpjm1    
     2263                        DO ji = 2, jpim1   ! NO vector opt.    
     2264                           zotx1(ji,jj) = 0.5  * ( un(ji,jj,1) + un(ji-1,jj,1) ) * zfr_l(ji,jj)   &       
     2265                                  &       + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &    
     2266                                  &                + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)    
     2267                           zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1) + vn(ji,jj-1,1) ) * zfr_l(ji,jj)   &    
     2268                                  &       + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &    
     2269                                  &                + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)    
     2270                        END DO    
     2271                     END DO    
     2272#if defined key_cice    
     2273                  ELSE    
     2274                     ! Temporarily Changed for HadGEM3    
     2275                     DO jj = 2, jpjm1    
     2276                        DO ji = 2, jpim1   ! NO vector opt.    
     2277                           zotx1(ji,jj) = (1.0-fr_iu(ji,jj)) * un(ji,jj,1)             &    
     2278                                &              + fr_iu(ji,jj) * 0.5 * ( u_ice(ji,jj-1) + u_ice(ji,jj) )     
     2279                           zoty1(ji,jj) = (1.0-fr_iv(ji,jj)) * vn(ji,jj,1)             &    
     2280                                &              + fr_iv(ji,jj) * 0.5 * ( v_ice(ji-1,jj) + v_ice(ji,jj) )     
     2281                        END DO    
     2282                     END DO    
     2283#endif    
     2284                  ENDIF  
    19422285               END SELECT 
    19432286               CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. ) 
    … …  
    19842327         IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components 
    19852328            !                                                                     ! Ocean component 
    1986             CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component  
    1987             CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component  
    1988             zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components  
    1989             zoty1(:,:) = ztmp2(:,:) 
    1990             IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component 
    1991                CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component  
    1992                CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component  
    1993                zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components  
    1994                zity1(:,:) = ztmp2(:,:) 
    1995             ENDIF 
     2329             IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN    
     2330                CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component    
     2331                CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component    
     2332                zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components    
     2333                zoty1(:,:) = ztmp2(:,:)    
     2334                IF( ssnd(jps_ivx1)%laction ) THEN                  ! Ice component    
     2335                   CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component    
     2336                   CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component    
     2337                   zitx1(:,:) = ztmp1(:,:) ! overwrite the components    
     2338                   zity1(:,:) = ztmp2(:,:)    
     2339                ENDIF    
     2340             ELSE    
     2341                ! Temporary code for HadGEM3 - will be removed eventually.    
     2342                ! Only applies when we want uvel on U grid and vvel on V grid    
     2343                ! Rotate U and V onto geographic grid before sending.    
     2344              
     2345                DO jj=2,jpjm1    
     2346                   DO ji=2,jpim1    
     2347                      ztmp1(ji,jj)=0.25*vmask(ji,jj,1)      &    
     2348                           *(zotx1(ji,jj)+zotx1(ji-1,jj)    &    
     2349                           +zotx1(ji,jj+1)+zotx1(ji-1,jj+1))    
     2350                      ztmp2(ji,jj)=0.25*umask(ji,jj,1)      &    
     2351                           *(zoty1(ji,jj)+zoty1(ji+1,jj)    &    
     2352                           +zoty1(ji,jj-1)+zoty1(ji+1,jj-1))    
     2353                   ENDDO    
     2354                ENDDO    
     2355                    
     2356                ! Ensure any N fold and wrap columns are updated    
     2357                CALL lbc_lnk(ztmp1, 'V', -1.0)    
     2358                CALL lbc_lnk(ztmp2, 'U', -1.0)    
     2359                    
     2360                ikchoix = -1    
     2361                CALL repcmo(zotx1,ztmp2,ztmp1,zoty1,zotx1,zoty1,ikchoix)    
     2362            ENDIF    
    19962363         ENDIF 
    19972364         ! 
    … …  
    20192386      ENDIF 
    20202387      ! 
     2388      !                                                      ! ------------------------- !    
     2389      !                                                      !  Surface current to waves !    
     2390      !                                                      ! ------------------------- !    
     2391      IF( ssnd(jps_ocxw)%laction .OR. ssnd(jps_ocyw)%laction ) THEN    
     2392          !        
     2393          !                                                  j+1  j     -----V---F    
     2394          ! surface velocity always sent from T point                    !       |    
     2395          !                                                       j      |   T   U    
     2396          !                                                              |       |    
     2397          !                                                   j   j-1   -I-------|    
     2398          !                                               (for I)        |       |    
     2399          !                                                             i-1  i   i    
     2400          !                                                              i      i+1 (for I)    
     2401          SELECT CASE( TRIM( sn_snd_crtw%cldes ) )    
     2402          CASE( 'oce only'             )      ! C-grid ==> T    
     2403             DO jj = 2, jpjm1    
     2404                DO ji = fs_2, fs_jpim1   ! vector opt.    
     2405                   zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )    
     2406                   zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji , jj-1,1) )     
     2407                END DO    
     2408             END DO    
     2409          CASE( 'weighted oce and ice' )       
     2410             SELECT CASE ( cp_ice_msh )    
     2411             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T    
     2412                DO jj = 2, jpjm1    
     2413                   DO ji = fs_2, fs_jpim1   ! vector opt.    
     2414                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un (ji-1,jj  ,1) ) * zfr_l(ji,jj)      
     2415                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn (ji  ,jj-1,1) ) * zfr_l(ji,jj)    
     2416                      zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)    
     2417                      zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)    
     2418                   END DO    
     2419                END DO    
     2420             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T    
     2421                DO jj = 2, jpjm1    
     2422                   DO ji = 2, jpim1   ! NO vector opt.    
     2423                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)      
     2424                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)      
     2425                      zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &    
     2426                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)    
     2427                      zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &    
     2428                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)    
     2429                   END DO    
     2430                END DO    
     2431             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T    
     2432                DO jj = 2, jpjm1    
     2433                   DO ji = 2, jpim1   ! NO vector opt.    
     2434                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)      
     2435                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)      
     2436                      zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &    
     2437                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)    
     2438                      zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &    
     2439                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)    
     2440                   END DO    
     2441                END DO    
     2442             END SELECT    
     2443             CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )    
     2444          CASE( 'mixed oce-ice'        )    
     2445             SELECT CASE ( cp_ice_msh )    
     2446             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T    
     2447                DO jj = 2, jpjm1    
     2448                   DO ji = fs_2, fs_jpim1   ! vector opt.    
     2449                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &    
     2450                         &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)    
     2451                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   &    
     2452                         &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)    
     2453                   END DO    
     2454                END DO    
     2455             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T    
     2456                DO jj = 2, jpjm1    
     2457                   DO ji = 2, jpim1   ! NO vector opt.    
     2458                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &       
     2459                         &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &    
     2460                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)    
     2461                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   &     
     2462                         &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &    
     2463                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)    
     2464                   END DO    
     2465                END DO    
     2466             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T    
     2467                DO jj = 2, jpjm1    
     2468                   DO ji = 2, jpim1   ! NO vector opt.    
     2469                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &       
     2470                         &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &    
     2471                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)    
     2472                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   &     
     2473                         &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &    
     2474                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)    
     2475                   END DO    
     2476                END DO    
     2477             END SELECT    
     2478          END SELECT    
     2479         CALL lbc_lnk( zotx1, ssnd(jps_ocxw)%clgrid, -1. )   ; CALL lbc_lnk( zoty1, ssnd(jps_ocyw)%clgrid, -1. )    
     2480         !    
     2481         !    
     2482         IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components    
     2483         !                                                                        ! Ocean component    
     2484            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->e', ztmp1 )       ! 1st component     
     2485            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->n', ztmp2 )       ! 2nd component     
     2486            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components     
     2487            zoty1(:,:) = ztmp2(:,:)     
     2488            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component    
     2489               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component     
     2490               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component     
     2491               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components     
     2492               zity1(:,:) = ztmp2(:,:)    
     2493            ENDIF    
     2494         ENDIF    
     2495         !    
     2496!         ! spherical coordinates to cartesian -> 2 components to 3 components    
     2497!         IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN    
     2498!            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents    
     2499!            ztmp2(:,:) = zoty1(:,:)    
     2500!            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )    
     2501!            !    
     2502!            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities    
     2503!               ztmp1(:,:) = zitx1(:,:)    
     2504!               ztmp1(:,:) = zity1(:,:)    
     2505!               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )    
     2506!            ENDIF    
     2507!         ENDIF    
     2508         !    
     2509         IF( ssnd(jps_ocxw)%laction )   CALL cpl_snd( jps_ocxw, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid    
     2510         IF( ssnd(jps_ocyw)%laction )   CALL cpl_snd( jps_ocyw, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid    
     2511         !     
     2512      ENDIF    
     2513      !    
     2514      IF( ssnd(jps_ficet)%laction ) THEN    
     2515         CALL cpl_snd( jps_ficet, isec, RESHAPE ( fr_i, (/jpi,jpj,1/) ), info )    
     2516      END IF    
     2517      !                                                      ! ------------------------- !    
     2518      !                                                      !   Water levels to waves   !    
     2519      !                                                      ! ------------------------- !    
     2520      IF( ssnd(jps_wlev)%laction ) THEN    
     2521         IF( ln_apr_dyn ) THEN     
     2522            IF( kt /= nit000 ) THEN     
     2523               ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )     
     2524            ELSE     
     2525               ztmp1(:,:) = sshb(:,:)     
     2526            ENDIF     
     2527         ELSE     
     2528            ztmp1(:,:) = sshn(:,:)     
     2529         ENDIF     
     2530         CALL cpl_snd( jps_wlev  , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )    
     2531      END IF 
    20212532      ! 
    20222533      !  Fields sent by OPA to SAS when doing OPA<->SAS coupling 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90

    r8918 r11286  
    2020   USE iom             ! IOM library 
    2121   USE in_out_manager  ! I/O manager 
     22   USE sbcwave         ! wave physics 
    2223   USE lib_mpp         ! distribued memory computing library 
    2324   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    … …  
    8788      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3 
    8889      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient 
     90      REAL(wp) ::   totwind               ! UKMO SHELF: Module of wind speed 
    8991      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables 
    9092      REAL     ::   cs                    ! UKMO SHELF: Friction co-efficient at surface 
    … …  
    115117         IF(lwm) WRITE ( numond, namsbc_flx )  
    116118         ! 
     119         IF(lwp) THEN                        ! Namelist print  
     120            WRITE(numout,*)  
     121            WRITE(numout,*) 'sbc_flx : Flux forcing'  
     122            WRITE(numout,*) '~~~~~~~~~~~'  
     123            WRITE(numout,*) '       Namelist namsbc_flx : shelf seas configuration (force with winds instead of momentum)'  
     124            WRITE(numout,*) '          shelf seas configuration    ln_shelf_flx    = ', ln_shelf_flx  
     125            WRITE(numout,*) '          relative wind speed         ln_rel_wind     = ', ln_rel_wind  
     126            WRITE(numout,*) '          wind multiplication factor  rn_wfac         = ', rn_wfac  
     127         ENDIF 
    117128         !                                         ! check: do we plan to use ln_dm2dc with non-daily forcing? 
    118129         IF( ln_dm2dc .AND. sn_qsr%nfreqh /= 24 )   & 
    … …  
    210221            END DO 
    211222         END DO 
     223         !                                                        ! add modification due to drag coefficient read from wave forcing  
     224         !                                                        ! this code is inefficient but put here to allow merging with another UKMO branch  
     225         IF( ln_shelf_flx ) THEN  
     226            ! calculate first the wind module, as it will be used later  
     227            DO jj = 2, jpjm1  
     228               DO ji = fs_2, fs_jpim1    ! vect. opt.  
     229                  ztx = zwnd_i(ji-1,jj  ) + zwnd_i(ji,jj)  
     230                  zty = zwnd_j(ji  ,jj-1) + zwnd_j(ji,jj)  
     231                  wndm(ji,jj) = 0.5 * SQRT( ztx * ztx + zty * zty )  
     232               END DO  
     233            END DO  
     234            CALL lbc_lnk( wndm(:,:), 'T', 1. )  
     235         
     236            IF( ln_cdgw .AND. nn_drag == jp_std ) THEN  
     237               IF( cpl_wdrag ) THEN   
     238                  ! reset utau and vtau to the wind components: the momentum will  
     239                  ! be calculated from the coupled value of the drag coefficient  
     240                  DO jj = 1, jpj  
     241                     DO ji = 1, jpi  
     242                        utau(ji,jj) = zwnd_i(ji,jj)  
     243                        vtau(ji,jj) = zwnd_j(ji,jj)  
     244                     END DO  
     245                  END DO  
     246               ELSE  
     247                  DO jj = 1, jpjm1  
     248                     DO ji = 1, jpim1  
     249                        utau(ji,jj) = zrhoa * 0.5 * ( cdn_wave(ji,jj) + cdn_wave(ji+1,jj) ) * zwnd_i(ji,jj) * &  
     250                                                                        0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) )  
     251                        vtau(ji,jj) = zrhoa * 0.5 * ( cdn_wave(ji,jj) + cdn_wave(ji,jj+1) ) * zwnd_j(ji,jj) * &  
     252                                                                        0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) )  
     253                     END DO  
     254                  END DO  
     255                  CALL lbc_lnk_multi( utau(:,:), 'U', -1., vtau(:,:), 'V', -1. )  
     256               ENDIF  
     257            ELSE IF( nn_drag == jp_const ) THEN  
     258               DO jj = 1, jpjm1  
     259                  DO ji = 1, jpim1  
     260                     utau(ji,jj) = zrhoa * zcdrag * zwnd_i(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) )  
     261                     vtau(ji,jj) = zrhoa * zcdrag * zwnd_j(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) )  
     262                  END DO  
     263               END DO  
     264               CALL lbc_lnk_multi( utau(:,:), 'U', -1., vtau(:,:), 'V', -1. )  
     265            ENDIF  
     266         ENDIF 
    212267         !                                                        ! add to qns the heat due to e-p 
    213268         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST 
    … …  
    230285               zmod = 0.5 * SQRT( ztx * ztx + zty * zty ) 
    231286               taum(ji,jj) = zmod 
     287               IF ( .NOT. (ln_shelf_flx .AND. ln_cpl)) THEN 
    232288               wndm(ji,jj) = SQRT( zmod * zcoef ) 
     289               ENDIF 
    233290            END DO 
    234291         END DO 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

    r8058 r11286  
    284284      CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 ) 
    285285      ! 
     286      ! In coupled mode get extra fields from CICE for passing back to atmosphere    
     287      IF ( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(nit000)    
     288      !     
    286289      IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_init') 
    287290      ! 
    … …  
    708711      IF( nn_timing == 1 )  CALL timing_start('cice_sbc_hadgam') 
    709712      ! 
    710       IF( kt == nit000 )  THEN 
    711          IF(lwp) WRITE(numout,*)'cice_sbc_hadgam' 
    712          IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' ) 
    713       ENDIF 
    714  
    715713      !                                         ! =========================== ! 
    716714      !                                         !   Prepare Coupling fields   ! 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r8058 r11286  
    8888      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl,   & 
    8989         &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   & 
    90          &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   & 
    91          &             nn_lsm    , nn_limflx , nn_components, ln_cpl 
     90         &             ln_ssr    , nn_isf    , nn_fwb, ln_wave, nn_lsm     , nn_limflx,   &  
     91                       nn_components, ln_cpl , ln_wavcpl, nn_drag 
    9292      INTEGER  ::   ios 
    9393      INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm 
    94       LOGICAL  ::   ll_purecpl 
    9594      !!---------------------------------------------------------------------- 
    9695 
    … …  
    131130         WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs  = ', ln_blk_mfs 
    132131         WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl      = ', ln_cpl 
    133          WRITE(numout,*) '              forced-coupled mixed formulation           ln_mixcpl   = ', ln_mixcpl 
     132         WRITE(numout,*) '              forced-coupled atm mixed formulation       ln_mixcpl   = ', ln_mixcpl  
     133         WRITE(numout,*) '              forced-coupled wav mixed formulation       ln_wavcpl   = ', ln_wavcpl 
     134         WRITE(numout,*) '              wave physics                               ln_wave     = ', ln_wave   
    134135         WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis    = ', lk_oasis 
    135136         WRITE(numout,*) '              components of your executable              nn_components = ', nn_components 
    136137         WRITE(numout,*) '              Multicategory heat flux formulation (LIM3) nn_limflx   = ', nn_limflx 
     138         WRITE(numout,*) '              momentum formulation                       nn_drag     = ', nn_drag 
    137139         WRITE(numout,*) '           Misc. options of sbc : ' 
    138140         WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn 
    … …  
    164166      IF ( nn_components == jp_iam_opa .AND. ln_cpl )   & 
    165167         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' ) 
    166       IF ( nn_components == jp_iam_opa .AND. ln_mixcpl )   & 
    167          &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' ) 
     168      IF ( nn_components == jp_iam_opa .AND. ( ln_mixcpl .OR. ln_wavcpl) )   &  
     169         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl or ln_wavcpl = T in OPA' ) 
    168170      IF ( ln_cpl .AND. .NOT. lk_oasis )    & 
    169171         &      CALL ctl_stop( 'STOP', 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' ) 
    170       IF( ln_mixcpl .AND. .NOT. lk_oasis )    & 
    171          &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires the cpp key key_oasis3' ) 
    172       IF( ln_mixcpl .AND. .NOT. ln_cpl )    & 
    173          &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires ln_cpl = T' ) 
    174       IF( ln_mixcpl .AND. nn_components /= jp_iam_nemo )    & 
    175          &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) is not yet working with sas-opa coupling via oasis' ) 
     172      IF( ( ln_mixcpl .OR. ln_wavcpl ) .AND. .NOT. lk_oasis )    &  
     173         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl or ln_wavcpl) requires the cpp key key_oasis3' )  
     174      IF( ( ln_mixcpl .OR. ln_wavcpl ) .AND. .NOT. ln_cpl )    &  
     175         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl or ln_wavcpl) requires ln_cpl = T' )  
     176      IF( ( ln_mixcpl .OR. ln_wavcpl ) .AND. nn_components /= jp_iam_nemo )    &  
     177         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl or ln_wavcpl) is not yet working with sas-opa coupling via oasis' ) 
    176178 
    177179      !                              ! allocate sbc arrays 
    … …  
    214216      IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk_core ) .AND. nn_components /= jp_iam_opa )   & 
    215217         &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' ) 
    216        
    217       IF ( ln_wave ) THEN 
    218       !Activated wave module but neither drag nor stokes drift activated 
    219          IF ( .NOT.(ln_cdgw .OR. ln_sdw) )   THEN 
    220             CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' ) 
    221       !drag coefficient read from wave model definable only with mfs bulk formulae and core  
    222          ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN        
    223              CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core') 
    224          ENDIF 
    225       ELSE 
    226       IF ( ln_cdgw .OR. ln_sdw  )                                         &  
    227          &   CALL ctl_stop('Not Activated Wave Module (ln_wave=F) but     & 
    228          & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ') 
    229       ENDIF  
    230218      !                          ! Choice of the Surface Boudary Condition (set nsbc) 
    231       ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl 
     219      ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl .AND. .NOT. ln_wavcpl 
    232220      ! 
    233221      icpt = 0 
    … …  
    261249         IF( nsbc == jp_mfs     )   WRITE(numout,*) '              MFS Bulk formulation' 
    262250         IF( nsbc == jp_none    )   WRITE(numout,*) '              OPA coupled to SAS via oasis' 
    263          IF( ln_mixcpl          )   WRITE(numout,*) '              + forced-coupled mixed formulation' 
     251         IF( ln_mixcpl          )   WRITE(numout,*) '              + forced-coupled mixed atm formulation'  
     252         IF( ln_wavcpl          )   WRITE(numout,*) '              + forced-coupled mixed wav formulation' 
    264253         IF( nn_components/= jp_iam_nemo )  & 
    265254            &                       WRITE(numout,*) '              + OASIS coupled SAS' 
    266255      ENDIF 
    267256      ! 
    268       IF( lk_oasis )   CALL sbc_cpl_init (nn_ice)   ! OASIS initialisation. must be done before: (1) first time step 
    269       !                                                     !                                            (2) the use of nn_fsbc 
    270  
     257      IF( lk_oasis ) THEN    
     258         IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )             
     259         CALL sbc_cpl_init (nn_ice)   ! OASIS initialisation. must be done before: (1) first time step    
     260                                      ! (2) the use of nn_fsbc    
     261      ENDIF    
    271262!     nn_fsbc initialization if OPA-SAS coupling via OASIS 
    272263!     sas model time step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly 
    … …  
    305296 
    306297      IF( nn_ice == 4      )   CALL cice_sbc_init( nsbc )      ! CICE initialisation 
     298      !   
     299      IF( ln_wave          )   CALL sbc_wave_init              ! surface wave initialisation   
     300      ! 
    307301       
    308302   END SUBROUTINE sbc_init 
    … …  
    325319      !!              - updte the ice fraction : fr_i 
    326320      !!---------------------------------------------------------------------- 
     321      USE bdydta, ONLY: bdy_dta   
     322      ! 
    327323      INTEGER, INTENT(in) ::   kt       ! ocean time step 
    328324      !!--------------------------------------------------------------------- 
    … …  
    345341      !                                            ! ---------------------------------------- ! 
    346342      ! 
    347       IF ( .NOT. lk_bdy ) then 
    348          IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc 
    349       ENDIF 
     343 
     344      IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc 
    350345                                                         ! (caution called before sbc_ssm) 
    351346      ! 
    … …  
    382377      END SELECT 
    383378 
    384       IF( ln_mixcpl )        CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing 
    385  
     379      IF( ln_mixcpl .OR. ln_wavcpl )  CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing  
     380       
     381      IF( ln_wave .AND. (ln_tauoc .OR. ln_tauw) ) CALL sbc_stress( )    ! Wave stress update   
     382      IF( lk_bdy )           CALL bdy_dta ( kt, time_offset=+1 )        ! update dynamic & tracer data at open boundaries  
     383                                                                        ! (caution called after sbc_ssm[_cpl] and before ice) 
    386384 
    387385      !                                            !==  Misc. Options  ==! 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r8058 r11286  
    8282      ALLOCATE( rnfmsk(jpi,jpj)         , rnfmsk_z(jpk)          ,     & 
    8383         &      h_rnf (jpi,jpj)         , nk_rnf  (jpi,jpj)      ,     & 
    84          &      rnf_tsc_b(jpi,jpj,jpts) , rnf_tsc (jpi,jpj,jpts) , STAT=sbc_rnf_alloc ) 
     84         &      rnf_tsc_b(jpi,jpj,jpts) , rnf_tsc (jpi,jpj,jpts) ,     &  
     85         &      xrnfmask(jpi,jpj,1)     , STAT=sbc_rnf_alloc ) 
    8586         ! 
    8687      IF( lk_mpp            )   CALL mpp_sum ( sbc_rnf_alloc ) 
    … …  
    128129      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
    129130         ! 
    130          IF( .NOT. l_rnfcpl )   rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )       ! updated runoff value at time step kt 
     131         IF( .NOT. l_rnfcpl )   &  
     132             rnf(:,:) = rnf(:,:) * (1. - xrnfmask(:,:,1)) + rn_rfact * sf_rnf(1)%fnow(:,:,1) * xrnfmask(:,:,1)  ! updated runoff value at time step kt 
    131133         ! 
    132134         !                                                     ! set temperature & salinity content of runoffs 
    … …  
    442444      ENDIF 
    443445      ! 
     446      xrnfmask(:,:,:) = 1.    ! default value for points using river forcing  
     447      IF (ln_usernfmask) THEN  
     448         IF(lwp) WRITE(numout,*)  
     449         IF(lwp) WRITE(numout,*) '          runoff mask read in a file'  
     450         CALL iom_open( 'rnfmask', inum )  
     451         CALL iom_get( inum, jpdom_data, 'rnfmask', xrnfmask(:,:,1), 1)  
     452         CALL iom_close( inum )  
     453      ENDIF  
     454      ! 
    444455      rnf(:,:) =  0._wp                         ! runoff initialisation 
    445456      rnf_tsc(:,:,:) = 0._wp                    ! runoffs temperature & salinty contents initilisation 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90

    r8058 r11286  
    2626 
    2727   PUBLIC   sbc_ssm         ! routine called by step.F90 
     28   PUBLIC   sbc_ssm_cpl     ! routine called by sbccpl.F90 
    2829   PUBLIC   sbc_ssm_init    ! routine called by sbcmod.F90 
    2930 
    … …  
    7778         sss_m(:,:) = zts(:,:,jp_sal) 
    7879         !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics) 
    79          IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
    80          ELSE                    ;   ssh_m(:,:) = sshn(:,:) 
     80         IF( .NOT. cpl_mslp ) THEN  
     81            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )  
     82            ELSE                    ;   ssh_m(:,:) = sshn(:,:)  
     83            ENDIF 
    8184         ENDIF 
    8285         ! 
    … …  
    99102            sss_m(:,:) = zcoef * zts(:,:,jp_sal) 
    100103            !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics) 
    101             IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) ) 
    102             ELSE                    ;   ssh_m(:,:) = zcoef * sshn(:,:) 
     104            IF( .NOT. cpl_mslp ) THEN  
     105               IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) )  
     106               ELSE                    ;   ssh_m(:,:) = zcoef * sshn(:,:)  
     107               ENDIF 
    103108            ENDIF 
    104109            ! 
    … …  
    113118            sst_m(:,:) = 0.e0 
    114119            sss_m(:,:) = 0.e0 
    115             ssh_m(:,:) = 0.e0 
     120            IF( .NOT. cpl_mslp ) ssh_m(:,:) = 0.e0 
    116121            IF( lk_vvl )   e3t_m(:,:) = 0.e0 
    117122            frq_m(:,:) = 0.e0 
    … …  
    127132         sss_m(:,:) = sss_m(:,:) + zts(:,:,jp_sal) 
    128133         !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics) 
    129          IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
    130          ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) 
     134         IF( .NOT. cpl_mslp ) THEN  
     135            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )  
     136            ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:)  
     137            ENDIF 
    131138         ENDIF 
    132139         ! 
    … …  
    143150            ssu_m(:,:) = ssu_m(:,:) * zcoef           ! mean suface current  [m/s] 
    144151            ssv_m(:,:) = ssv_m(:,:) * zcoef           ! 
    145             ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH             [m] 
     152            IF( .NOT. cpl_mslp ) ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH             [m] 
    146153            IF( lk_vvl )   e3t_m(:,:) = fse3t_m(:,:) * zcoef   ! mean vertical scale factor [m] 
    147154            frq_m(:,:) = frq_m(:,:) * zcoef   ! mean fraction of solar net radiation absorbed in the 1st T level [-] 
    … …  
    161168            CALL iom_rstput( kt, nitrst, numrow, 'sst_m'  , sst_m  ) 
    162169            CALL iom_rstput( kt, nitrst, numrow, 'sss_m'  , sss_m  ) 
    163             CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  ) 
     170            IF( .NOT. cpl_mslp ) CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  ) 
    164171            IF( lk_vvl )   CALL iom_rstput( kt, nitrst, numrow, 'e3t_m'  , e3t_m  ) 
    165172            CALL iom_rstput( kt, nitrst, numrow, 'frq_m'  , frq_m  ) 
    … …  
    174181         CALL iom_put( 'sst_m', sst_m ) 
    175182         CALL iom_put( 'sss_m', sss_m ) 
    176          CALL iom_put( 'ssh_m', ssh_m ) 
     183         IF( .NOT. cpl_mslp ) CALL iom_put( 'ssh_m', ssh_m ) 
    177184         IF( lk_vvl )   CALL iom_put( 'e3t_m', e3t_m ) 
    178185         CALL iom_put( 'frq_m', frq_m ) 
    … …  
    180187      ! 
    181188   END SUBROUTINE sbc_ssm 
     189 
     190   SUBROUTINE sbc_ssm_cpl( kt )  
     191      !!---------------------------------------------------------------------  
     192      !!                   ***  ROUTINE sbc_ssm_cpl  ***  
     193      !!                       
     194      !! ** Purpose :   provide ocean surface variable to sea-surface boundary  
     195      !!                condition computation when pressure is read from coupling  
     196      !!                  
     197      !! ** Method  :   The inverse barometer ssh (i.e. ssh associated with Patm)  
     198      !!                is added to ssh_m when ln_apr_dyn = T. Required for sea-ice dynamics.  
     199      !!---------------------------------------------------------------------  
     200      INTEGER, INTENT(in) ::   kt   ! ocean time step  
     201      !  
     202      REAL(wp) ::   zcoef       ! local scalar  
     203      !!---------------------------------------------------------------------  
     204      !  
     205      IF( nn_fsbc == 1 ) THEN                             !      Instantaneous surface fields        !  
     206         !                                                ! ---------------------------------------- !  
     207         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )  
     208         ELSE                    ;   ssh_m(:,:) = sshn(:,:)  
     209         ENDIF  
     210      ELSE  
     211         !                                                ! ----------------------------------------------- !  
     212         IF( kt == nit000 .AND. .NOT. l_ssm_mean ) THEN   !   Initialisation: 1st time-step, no input means !  
     213            !                                             ! ----------------------------------------------- !  
     214            IF(lwp) WRITE(numout,*)  
     215            IF(lwp) WRITE(numout,*) '~~~~~~~   mean ssh field initialised to instantaneous values'  
     216            zcoef = REAL( nn_fsbc - 1, wp )  
     217            zcoef = REAL( nn_fsbc - 1, wp )  
     218            IF( ln_apr_dyn ) THEN    ;  ssh_m(:,:) = zcoef * ( sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) )  
     219            ELSE                     ;  ssh_m(:,:) = zcoef * sshn(:,:)  
     220            ENDIF  
     221            !                                             ! ---------------------------------------- !  
     222         ELSEIF( MOD( kt - 2 , nn_fsbc ) == 0 ) THEN      !   Initialisation: New mean computation   !  
     223            !                                             ! ---------------------------------------- !  
     224            ssh_m(:,:) = 0.e0  
     225         ENDIF  
     226    
     227         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )  
     228         ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:)  
     229         ENDIF  
     230         !                                                ! ---------------------------------------- !  
     231         IF( MOD( kt - 1 , nn_fsbc ) == 0 ) THEN          !   Mean value at each nn_fsbc time-step   !  
     232            !                                             ! ---------------------------------------- !  
     233            zcoef = 1. / REAL( nn_fsbc, wp )  
     234            ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH [m]  
     235         ENDIF  
     236         !                                                ! ---------------------------------------- !  
     237         IF( lrst_oce ) THEN                              !      Write in the ocean restart file     !  
     238            !                                             ! ---------------------------------------- !  
     239            IF(lwp) WRITE(numout,*)  
     240            IF(lwp) WRITE(numout,*) 'sbc_ssm_cpl : ssh mean field written in ocean restart file ',   &  
     241               &                    'at it= ', kt,' date= ', ndastp  
     242            IF(lwp) WRITE(numout,*) '~~~~~~~'  
     243            CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  )  
     244         ENDIF  
     245      ENDIF  
     246      !  
     247      IF( MOD( kt - 1 , nn_fsbc ) == 0 ) THEN          !   Mean value at each nn_fsbc time-step   !  
     248         CALL iom_put( 'ssh_m', ssh_m )  
     249      ENDIF  
     250      !  
     251   END SUBROUTINE sbc_ssm_cpl 
    182252 
    183253   SUBROUTINE sbc_ssm_init 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r8058 r11286  
    44   !! Wave module  
    55   !!====================================================================== 
    6    !! History :  3.3.1  !   2011-09  (Adani M)  Original code: Drag Coefficient  
    7    !!         :  3.4    !   2012-10  (Adani M)                 Stokes Drift  
     6   !! History :  3.3  !   2011-09  (Adani M)  Original code: Drag Coefficient  
     7   !!         :  3.4  !   2012-10  (Adani M)                 Stokes Drift  
     8   !!            3.6  !   2014-09  (Clementi E, Oddo P)New Stokes Drift Computation 
     9   !!             -   !   2016-12  (G. Madec, E. Clementi) update Stoke drift computation  
     10   !!                                                     + add sbc_wave_ini routine 
    811   !!---------------------------------------------------------------------- 
    9    USE iom             ! I/O manager library 
    10    USE in_out_manager  ! I/O manager 
    11    USE lib_mpp         ! distribued memory computing library 
    12    USE fldread        ! read input fields 
    13    USE oce 
    14    USE sbc_oce        ! Surface boundary condition: ocean fields 
    15    USE domvvl 
    16  
    17     
     12 
    1813   !!---------------------------------------------------------------------- 
    19    !!   sbc_wave       : read drag coefficient from wave model in netcdf files  
     14   !!   sbc_stokes    : calculate 3D Stokes-drift velocities 
     15   !!   sbc_wave      : wave data from wave model in netcdf files  
     16   !!   sbc_wave_init : initialisation fo surface waves 
    2017   !!---------------------------------------------------------------------- 
     18   USE oce            ! ocean variables 
     19   USE sbc_oce        ! Surface boundary condition: ocean fields 
     20   USE bdy_oce        ! open boundary condition variables 
     21   USE domvvl         ! domain: variable volume layers 
     22   ! 
     23   USE iom            ! I/O manager library 
     24   USE in_out_manager ! I/O manager 
     25   USE lib_mpp        ! distribued memory computing library 
     26   USE fldread        ! read input fields 
     27   USE wrk_nemo       ! 
     28   USE phycst         ! physical constants  
    2129 
    2230   IMPLICIT NONE 
    2331   PRIVATE 
    2432 
    25    PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs 
     33   PUBLIC   sbc_stokes      ! routine called in sbccpl  
     34   PUBLIC   sbc_stress      ! routine called in sbcmod 
     35   PUBLIC   sbc_wave        ! routine called in sbcmod  
     36   PUBLIC   sbc_wave_init   ! routine called in sbcmod 
    2637    
    27    INTEGER , PARAMETER ::   jpfld  = 3           ! maximum number of files to read for srokes drift 
    28    INTEGER , PARAMETER ::   jp_usd = 1           ! index of stokes drift  (i-component) (m/s)    at T-point 
    29    INTEGER , PARAMETER ::   jp_vsd = 2           ! index of stokes drift  (j-component) (m/s)    at T-point 
    30    INTEGER , PARAMETER ::   jp_wn  = 3           ! index of wave number                 (1/m)    at T-point 
    31    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
    32    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
    33    REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave  
    34    REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum  
    35    REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d  
    36  
    37    !! * Substitutions 
    38 #  include "domzgr_substitute.h90" 
     38   ! Variables checking if the wave parameters are coupled (if not, they are read from file) 
     39   LOGICAL, PUBLIC ::   cpl_hsig   = .FALSE. 
     40   LOGICAL, PUBLIC ::   cpl_phioc  = .FALSE. 
     41   LOGICAL, PUBLIC ::   cpl_sdrft  = .FALSE. 
     42   LOGICAL, PUBLIC ::   cpl_wper   = .FALSE. 
     43   LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE. 
     44   LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE. 
     45   LOGICAL, PUBLIC ::   cpl_tauoc  = .FALSE. 
     46   LOGICAL, PUBLIC ::   cpl_tauw   = .FALSE. 
     47   LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE. 
     48 
     49   INTEGER ::   nn_sdrift      ! type of parameterization to calculate vertical Stokes drift 
     50   INTEGER, PARAMETER ::   jp_breivik  = 0     ! Breivik 2015: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] 
     51   INTEGER, PARAMETER ::   jp_phillips = 1     ! Phillips:     v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))] 
     52   INTEGER, PARAMETER ::   jp_peakph   = 2     ! Phillips using the peak wave number read from wave model instead of the inverse depth scale 
     53 
     54   INTEGER ::   jpfld    ! number of files to read for stokes drift 
     55   INTEGER ::   jp_usd   ! index of stokes drift  (i-component) (m/s)    at T-point 
     56   INTEGER ::   jp_vsd   ! index of stokes drift  (j-component) (m/s)    at T-point 
     57   INTEGER ::   jp_hsw   ! index of significant wave hight      (m)      at T-point 
     58   INTEGER ::   jp_wmp   ! index of mean wave period            (s)      at T-point 
     59   INTEGER ::   jp_wfr   ! index of wave peak frequency         (s^-1)   at T-point 
     60 
     61   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
     62   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
     63   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao 
     64   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 
     65   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_tauw ! structure of input fields (file informations, fields read) ocean stress components from wave model 
     66   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_phioc ! structure of input fields (file informations, fields read) wave to ocean energy  
     67   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !: 
     68   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw, wmp, wnum      !: 
     69   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   wfreq               !:  
     70   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   rn_crban            !: Craig and Banner constant for surface breaking waves mixing 
     71   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !: 
     72   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauw_x              !: 
     73   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauw_y              !: 
     74   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !: 
     75   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   div_sd              !: barotropic stokes drift divergence 
     76   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   ut0sd, vt0sd        !: surface Stokes drift velocities at t-point 
     77   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd  , vsd  , wsd   !: Stokes drift velocities at u-, v- & w-points, resp. 
     78 
     79#  include "vectopt_loop_substitute.h90" 
    3980   !!---------------------------------------------------------------------- 
    4081   !! NEMO/OPA 4.0 , NEMO Consortium (2011)  
    … …  
    4485CONTAINS 
    4586 
     87   SUBROUTINE sbc_stokes( ) 
     88      !!--------------------------------------------------------------------- 
     89      !!                     ***  ROUTINE sbc_stokes  *** 
     90      !! 
     91      !! ** Purpose :   compute the 3d Stokes Drift according to Breivik et al., 
     92      !!                2014 (DOI: 10.1175/JPO-D-14-0020.1) 
     93      !! 
     94      !! ** Method  : - Calculate Stokes transport speed  
     95      !!              - Calculate horizontal divergence  
     96      !!              - Integrate the horizontal divergenze from the bottom  
     97      !! ** action   
     98      !!--------------------------------------------------------------------- 
     99      INTEGER  ::   jj, ji, jk   ! dummy loop argument 
     100      INTEGER  ::   ik           ! local integer  
     101      REAL(wp) ::  ztransp, zfac, ztemp 
     102      REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v, zda_u, zda_v 
     103      REAL(wp), DIMENSION(:,:)  , POINTER ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd   ! 2D workspace 
     104      REAL(wp), DIMENSION(:,:,:), POINTER ::   ze3divh                            ! 3D workspace 
     105 
     106      !!--------------------------------------------------------------------- 
     107      ! 
     108 
     109      CALL wrk_alloc( jpi,jpj,jpk,   ze3divh ) 
     110      CALL wrk_alloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd )  
     111      ! 
     112      ! select parameterization for the calculation of vertical Stokes drift 
     113      ! exp. wave number at t-point 
     114      IF( nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips ) THEN   ! (Eq. (19) in Breivick et al. (2014) ) 
     115         zfac = 2.0_wp * rpi / 16.0_wp 
     116         DO jj = 1, jpj                
     117            DO ji = 1, jpi 
     118               ! Stokes drift velocity estimated from Hs and Tmean 
     119               ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp ) 
     120               ! Stokes surface speed 
     121               tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj)) 
     122               ! Wavenumber scale 
     123               zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp ) 
     124            END DO 
     125         END DO 
     126         DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points 
     127            DO ji = 1, jpim1 
     128               zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) 
     129               zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) 
     130               ! 
     131               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 
     132               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 
     133            END DO 
     134         END DO 
     135      ELSE IF( nn_sdrift==jp_peakph ) THEN    ! peak wave number calculated from the peak frequency received by the wave model 
     136         DO jj = 1, jpjm1               
     137            DO ji = 1, jpim1 
     138               zk_u(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji+1,jj)*wfreq(ji+1,jj) ) / grav 
     139               zk_v(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji,jj+1)*wfreq(ji,jj+1) ) / grav 
     140               ! 
     141               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 
     142               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 
     143            END DO 
     144         END DO 
     145      ENDIF 
     146      ! 
     147      !                       !==  horizontal Stokes Drift 3D velocity  ==! 
     148      IF( nn_sdrift==jp_breivik ) THEN 
     149         DO jk = 1, jpkm1 
     150            DO jj = 2, jpjm1 
     151               DO ji = 2, jpim1 
     152                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 
     153                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 
     154                  !                           
     155                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth 
     156                  zkh_v = zk_v(ji,jj) * zdep_v 
     157                  !                                ! Depth attenuation 
     158                  zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u ) 
     159                  zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v ) 
     160                  ! 
     161                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 
     162                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 
     163               END DO 
     164            END DO 
     165         END DO 
     166      ELSE IF( nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph ) THEN 
     167         DO jk = 1, jpkm1 
     168            DO jj = 2, jpjm1 
     169               DO ji = 2, jpim1 
     170                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 
     171                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 
     172                  !                           
     173                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth 
     174                  zkh_v = zk_v(ji,jj) * zdep_v 
     175                  !                                ! Depth attenuation 
     176                  zda_u = EXP( -2.0_wp*zkh_u ) - SQRT(2.0_wp*rpi*zkh_u) * ERFC(SQRT(2.0_wp*zkh_u)) 
     177                  zda_v = EXP( -2.0_wp*zkh_v ) - SQRT(2.0_wp*rpi*zkh_v) * ERFC(SQRT(2.0_wp*zkh_v)) 
     178                  ! 
     179                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 
     180                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 
     181               END DO 
     182            END DO 
     183         END DO 
     184      ENDIF 
     185 
     186      CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. ) 
     187      ! 
     188      !                       !==  vertical Stokes Drift 3D velocity  ==! 
     189      ! 
     190      DO jk = 1, jpkm1               ! Horizontal e3*divergence 
     191         DO jj = 2, jpj 
     192            DO ji = fs_2, jpi 
     193               ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * usd(ji,  jj,jk)    & 
     194                  &                 - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * usd(ji-1,jj,jk)    & 
     195                  &                 + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vsd(ji,jj  ,jk)    & 
     196                  &                 - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vsd(ji,jj-1,jk)  ) * r1_e12t(ji,jj) 
     197            END DO 
     198         END DO 
     199      END DO 
     200      ! 
     201      IF( .NOT. AGRIF_Root() ) THEN 
     202         IF( nbondi ==  1 .OR. nbondi == 2 )   ze3divh(nlci-1,   :  ,:) = 0._wp      ! east 
     203         IF( nbondi == -1 .OR. nbondi == 2 )   ze3divh(  2   ,   :  ,:) = 0._wp      ! west 
     204         IF( nbondj ==  1 .OR. nbondj == 2 )   ze3divh(  :   ,nlcj-1,:) = 0._wp      ! north 
     205         IF( nbondj == -1 .OR. nbondj == 2 )   ze3divh(  :   ,  2   ,:) = 0._wp      ! south 
     206      ENDIF 
     207      ! 
     208      CALL lbc_lnk( ze3divh, 'T', 1. ) 
     209      ! 
     210      IF( .NOT. lk_vvl ) THEN   ;   ik = 1   ! none zero velocity through the sea surface 
     211      ELSE                      ;   ik = 2   ! w=0 at the surface (set one for all in sbc_wave_init) 
     212      ENDIF 
     213      DO jk = jpkm1, ik, -1          ! integrate from the bottom the hor. divergence (NB: at k=jpk w is always zero) 
     214         wsd(:,:,jk) = wsd(:,:,jk+1) - ze3divh(:,:,jk) 
     215      END DO 
     216#if defined key_bdy 
     217      IF( lk_bdy ) THEN 
     218         DO jk = 1, jpkm1 
     219            wsd(:,:,jk) = wsd(:,:,jk) * bdytmask(:,:) 
     220         END DO 
     221      ENDIF 
     222#endif 
     223      !                       !==  Horizontal divergence of barotropic Stokes transport  ==! 
     224      div_sd(:,:) = 0._wp 
     225      DO jk = 1, jpkm1                                 !  
     226        div_sd(:,:) = div_sd(:,:) + ze3divh(:,:,jk) 
     227      END DO 
     228      ! 
     229      CALL iom_put( "ustokes",  usd  ) 
     230      CALL iom_put( "vstokes",  vsd  ) 
     231      CALL iom_put( "wstokes",  wsd  ) 
     232      ! 
     233      CALL wrk_dealloc( jpi,jpj,jpk,   ze3divh ) 
     234      CALL wrk_dealloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 
     235      ! 
     236   END SUBROUTINE sbc_stokes 
     237 
     238 
     239   SUBROUTINE sbc_stress( ) 
     240      !!--------------------------------------------------------------------- 
     241      !!                     ***  ROUTINE sbc_stress  *** 
     242      !! 
     243      !! ** Purpose :   Updates the ocean momentum modified by waves 
     244      !! 
     245      !! ** Method  : - Calculate u,v components of stress depending on stress 
     246      !!                model  
     247      !!              - Calculate the stress module 
     248      !!              - The wind module is not modified by waves  
     249      !! ** action   
     250      !!--------------------------------------------------------------------- 
     251      INTEGER  ::   jj, ji   ! dummy loop argument 
     252      ! 
     253      IF( ln_tauoc ) THEN 
     254         utau(:,:) = utau(:,:)*tauoc_wave(:,:) 
     255         vtau(:,:) = vtau(:,:)*tauoc_wave(:,:) 
     256         taum(:,:) = taum(:,:)*tauoc_wave(:,:) 
     257      ENDIF 
     258      ! 
     259      IF( ln_tauw ) THEN 
     260         DO jj = 1, jpjm1 
     261            DO ji = 1, jpim1 
     262               ! Stress components at u- & v-points 
     263               utau(ji,jj) = 0.5_wp * ( tauw_x(ji,jj) + tauw_x(ji+1,jj) ) 
     264               vtau(ji,jj) = 0.5_wp * ( tauw_y(ji,jj) + tauw_y(ji,jj+1) ) 
     265               ! 
     266               ! Stress module at t points 
     267               taum(ji,jj) = SQRT( tauw_x(ji,jj)*tauw_x(ji,jj) + tauw_y(ji,jj)*tauw_y(ji,jj) ) 
     268            END DO 
     269         END DO 
     270         CALL lbc_lnk_multi( utau(:,:), 'U', -1. , vtau(:,:), 'V', -1. , taum(:,: ), 'T', -1. ) 
     271      ENDIF 
     272      ! 
     273   END SUBROUTINE sbc_stress 
     274 
     275 
    46276   SUBROUTINE sbc_wave( kt ) 
    47277      !!--------------------------------------------------------------------- 
    48       !!                     ***  ROUTINE sbc_apr  *** 
    49       !! 
    50       !! ** Purpose :   read drag coefficient from wave model  in netcdf files. 
     278      !!                     ***  ROUTINE sbc_wave  *** 
     279      !! 
     280      !! ** Purpose :   read wave parameters from wave model  in netcdf files. 
    51281      !! 
    52282      !! ** Method  : - Read namelist namsbc_wave 
    53283      !!              - Read Cd_n10 fields in netcdf files  
    54284      !!              - Read stokes drift 2d in netcdf files  
    55       !!              - Read wave number      in netcdf files  
    56       !!              - Compute 3d stokes drift using monochromatic 
    57       !! ** action  :    
    58       !!                
    59       !!--------------------------------------------------------------------- 
    60       USE oce,  ONLY : un,vn,hdivn,rotn 
    61       USE divcur 
    62       USE wrk_nemo 
    63 #if defined key_bdy 
    64       USE bdy_oce, ONLY : bdytmask 
    65 #endif 
    66       INTEGER, INTENT( in  ) ::  kt       ! ocean time step 
    67       INTEGER                ::  ierror   ! return error code 
    68       INTEGER                ::  ifpr, jj,ji,jk  
    69       INTEGER                ::   ios     ! Local integer output status for namelist read 
    70       REAL(wp),DIMENSION(:,:,:),POINTER             ::  udummy,vdummy,hdivdummy,rotdummy 
    71       REAL                                          ::  z2dt,z1_2dt 
    72       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read 
     285      !!              - Read wave number in netcdf files  
     286      !!              - Compute 3d stokes drift using Breivik et al.,2014 
     287      !!                formulation 
     288      !! ** action   
     289      !!--------------------------------------------------------------------- 
     290      INTEGER, INTENT(in   ) ::   kt   ! ocean time step 
     291      !!--------------------------------------------------------------------- 
     292      ! 
     293      IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN     !==  Neutral drag coefficient  ==! 
     294         CALL fld_read( kt, nn_fsbc, sf_cd )             ! read from external forcing 
     295         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
     296         ! check that the drag coefficient contains proper information even if 
     297         ! the masks do not match - the momentum stress is not masked! 
     298         WHERE( cdn_wave < 0.0 ) cdn_wave = 1.5e-3 
     299         WHERE( cdn_wave > 1.0 ) cdn_wave = 1.5e-3 
     300      ENDIF 
     301 
     302      IF( ln_tauoc .AND. .NOT. cpl_tauoc ) THEN    !==  Wave induced stress  ==! 
     303         CALL fld_read( kt, nn_fsbc, sf_tauoc )          ! read wave norm stress from external forcing 
     304         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 
     305         WHERE( tauoc_wave <   0.0 ) tauoc_wave = 1.0 
     306         WHERE( tauoc_wave > 100.0 ) tauoc_wave = 1.0 
     307      ENDIF 
     308 
     309      IF( ln_tauw .AND. .NOT. cpl_tauw ) THEN      !==  Wave induced stress  ==! 
     310         CALL fld_read( kt, nn_fsbc, sf_tauw )           ! read ocean stress components from external forcing (T grid) 
     311         tauw_x(:,:) = sf_tauw(1)%fnow(:,:,1) 
     312         WHERE( tauw_x < -100.0 ) tauw_x = 0.0 
     313         WHERE( tauw_x >  100.0 ) tauw_x = 0.0 
     314 
     315         tauw_y(:,:) = sf_tauw(2)%fnow(:,:,1) 
     316         WHERE( tauw_y < -100.0 ) tauw_y = 0.0 
     317         WHERE( tauw_y >  100.0 ) tauw_y = 0.0 
     318      ENDIF 
     319 
     320      IF( ln_phioc .AND. .NOT. cpl_phioc ) THEN    !==  Wave to ocean energy  ==! 
     321         CALL fld_read( kt, nn_fsbc, sf_phioc )          ! read wave to ocean energy from external forcing 
     322         rn_crban(:,:) = 29.0 * sf_phioc(1)%fnow(:,:,1)     ! ! Alfa is phioc*sqrt(rau0/zrhoa)  : rau0=water density, zhroa= air density 
     323         WHERE( rn_crban > 1.e8   ) rn_crban = 0.0    !remove first mask mistmatch points, then cap values in case of low friction velocity 
     324         WHERE( rn_crban < 0.0    ) rn_crban = 0.0 
     325         WHERE( rn_crban > 1000.0 ) rn_crban = 1000.0 
     326      ENDIF 
     327 
     328      IF( ln_sdw .OR. ln_rough )  THEN             !==  Computation of the 3d Stokes Drift  ==!  
     329         ! 
     330         IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled 
     331            CALL fld_read( kt, nn_fsbc, sf_sd )          ! read wave parameters from external forcing 
     332            IF( jp_hsw > 0 ) THEN 
     333               hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height 
     334               WHERE( hsw > 100.0 ) hsw = 0.0 
     335               WHERE( hsw <   0.0 ) hsw = 0.0 
     336            ENDIF 
     337            IF( jp_wmp > 0 ) THEN 
     338               wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
     339               WHERE( wmp > 100.0 ) wmp = 0.0 
     340               WHERE( wmp <   0.0 ) wmp = 0.0 
     341            ENDIF 
     342            IF( jp_wfr > 0 ) THEN 
     343               wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1)   ! Peak wave frequency  
     344               WHERE( wfreq <    0.0 ) wfreq = 0.0  
     345               WHERE( wfreq >  100.0 ) wfreq = 0.0 
     346            ENDIF 
     347            IF( jp_usd > 0 ) THEN 
     348               ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
     349               WHERE( ut0sd < -100.0 ) ut0sd = 1.0 
     350               WHERE( ut0sd >  100.0 ) ut0sd = 1.0 
     351            ENDIF 
     352            IF( jp_vsd > 0 ) THEN 
     353               vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
     354               WHERE( vt0sd < -100.0 ) vt0sd = 1.0 
     355               WHERE( vt0sd >  100.0 ) vt0sd = 1.0 
     356            ENDIF 
     357         ENDIF 
     358      ENDIF 
     359      ! 
     360      IF( ln_sdw ) THEN 
     361         ! Read also wave number if needed, so that it is available in coupling routines 
     362         IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN 
     363            CALL fld_read( kt, nn_fsbc, sf_wn )          ! read wave parameters from external forcing 
     364            wnum(:,:) = sf_wn(1)%fnow(:,:,1) 
     365         ENDIF 
     366            
     367         !                                         !==  Computation of the 3d Stokes Drift  ==!  
     368         ! 
     369         IF( ((nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) .AND. & 
     370                          jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0) .OR. & 
     371             (nn_sdrift==jp_peakph .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0) ) & 
     372            CALL sbc_stokes()            ! Calculate only if required fields are read 
     373         !                               ! In coupled wave model-NEMO case the call is done after coupling 
     374         ! 
     375      ENDIF 
     376      ! 
     377   END SUBROUTINE sbc_wave 
     378 
     379 
     380   SUBROUTINE sbc_wave_init 
     381      !!--------------------------------------------------------------------- 
     382      !!                     ***  ROUTINE sbc_wave_init  *** 
     383      !! 
     384      !! ** Purpose :   read wave parameters from wave model  in netcdf files. 
     385      !! 
     386      !! ** Method  : - Read namelist namsbc_wave 
     387      !!              - Read Cd_n10 fields in netcdf files  
     388      !!              - Read stokes drift 2d in netcdf files  
     389      !!              - Read wave number in netcdf files  
     390      !!              - Compute 3d stokes drift using Breivik et al.,2014 
     391      !!                formulation 
     392      !! ** action   
     393      !!--------------------------------------------------------------------- 
     394      INTEGER ::   ierror, ios   ! local integer 
     395      INTEGER ::   ifpr 
     396      !! 
    73397      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files 
    74       TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_wn   ! informations about the fields to be read 
    75       !!--------------------------------------------------------------------- 
    76       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn 
    77       !!--------------------------------------------------------------------- 
    78  
    79       !!---------------------------------------------------------------------- 
    80       ! 
    81       ! 
    82       !                                         ! -------------------- ! 
    83       IF( kt == nit000 ) THEN                   ! First call kt=nit000 ! 
    84          !                                      ! -------------------- ! 
    85          REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model 
    86          READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 
    87 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 
    88  
    89          REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 
    90          READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 
    91 902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp ) 
    92          IF(lwm) WRITE ( numond, namsbc_wave ) 
    93          ! 
    94  
    95          IF ( ln_cdgw ) THEN 
     398      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i, slf_j     ! array of namelist informations on the fields to read 
     399      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_phioc, & 
     400                             &   sn_hsw, sn_wmp, sn_wfr, sn_wnum , & 
     401                             &   sn_tauoc, sn_tauwx, sn_tauwy      ! information about the fields to be read 
     402      ! 
     403      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, sn_wnum, sn_phioc,          & 
     404                             sn_tauoc, sn_tauwx, sn_tauwy,                                                       & 
     405                             ln_cdgw, ln_sdw, ln_stcor, ln_phioc, ln_tauoc, ln_tauw, ln_zdfqiao, ln_rough,       & 
     406                             nn_sdrift, nn_wmix 
     407      !!--------------------------------------------------------------------- 
     408      ! 
     409      REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model 
     410      READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 
     411901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 
     412          
     413      REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 
     414      READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 
     415902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp ) 
     416      IF(lwm) WRITE ( numond, namsbc_wave ) 
     417      ! 
     418      IF(lwp) THEN               !* Parameter print 
     419         WRITE(numout,*) 
     420         WRITE(numout,*) 'sbc_wave_init: wave physics' 
     421         WRITE(numout,*) '~~~~~~~~' 
     422         WRITE(numout,*) '   Namelist namsbc_wave : set wave physics parameters' 
     423         WRITE(numout,*) '      Stokes drift corr. to vert. velocity    ln_sdw      = ', ln_sdw  
     424         WRITE(numout,*) '        vertical parametrization              nn_sdrift   = ', nn_sdrift  
     425         WRITE(numout,*) '      Stokes coriolis term                    ln_stcor    = ', ln_stcor  
     426         WRITE(numout,*) '      wave modified ocean stress              ln_tauoc    = ', ln_tauoc  
     427         WRITE(numout,*) '      wave modified ocean stress components   ln_tauw     = ', ln_tauw  
     428         WRITE(numout,*) '      wave to ocean energy                    ln_phioc    = ', ln_phioc 
     429         WRITE(numout,*) '        vertical mixing parametrization       nn_wmix     = ', nn_wmix  
     430         WRITE(numout,*) '      neutral drag coefficient                ln_cdgw     = ', ln_cdgw 
     431         WRITE(numout,*) '      wave roughness length modification      ln_rough    = ', ln_rough  
     432         WRITE(numout,*) '      Qiao vertical mixing formulation        ln_zdfqiao  = ', ln_zdfqiao 
     433      ENDIF 
     434 
     435      IF ( ln_wave ) THEN 
     436         ! Activated wave physics but no wave physics components activated  
     437         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_tauw .OR. ln_stcor .OR. ln_phioc & 
     438                                                                    .OR. ln_rough .OR. ln_zdfqiao) )   THEN   
     439             CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_tauw=F, ln_stcor=F ', & 
     440                                                      'ln_phioc=F, ln_rough=F, ln_zdfqiao=F' )  
     441         ELSE 
     442         IF (ln_stcor .AND. .NOT. ln_sdw) &  
     443             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 
     444         IF ( ln_cdgw .AND. .NOT.(nn_drag==jp_ukmo .OR. nn_drag==jp_std .OR. nn_drag==jp_const .OR. nn_drag==jp_mcore) ) &  
     445             CALL ctl_stop( 'The chosen nn_drag for momentum calculation must be 0, 1, 2, or 3') 
     446         IF ( ln_cdgw .AND. ln_blk_core .AND. nn_drag==0 ) & 
     447             CALL ctl_stop( 'The chosen nn_drag for momentum calculation in core forcing must be 1, 2, or 3') 
     448         IF ( ln_cdgw .AND. ln_flx .AND. nn_drag==3 ) & 
     449             CALL ctl_stop( 'The chosen nn_drag for momentum calculation in direct forcing must be 0, 1, or 2') 
     450         IF( ln_phioc .AND. .NOT.(nn_wmix==jp_craigbanner .OR. nn_wmix==jp_janssen) ) &  
     451            CALL ctl_stop( 'The chosen nn_wmix for wave vertical mixing must be 0, or 1' ) 
     452         IF( ln_sdw .AND. .NOT.(nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph) ) &  
     453            CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' ) 
     454         IF( ln_zdfqiao .AND. .NOT.ln_sdw ) &  
     455            CALL ctl_stop( 'Qiao vertical mixing can not be used without Stokes drift (ln_sdw)' ) 
     456         IF( ln_tauoc .AND. ln_tauw ) &  
     457            CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', & 
     458                                     '(ln_tauoc=.true. and ln_tauw=.true.)' ) 
     459         IF( ln_tauoc ) & 
     460             CALL ctl_warn( 'You are subtracting the wave stress to the ocean (ln_tauoc=.true.)' ) 
     461         IF( ln_tauw ) & 
     462             CALL ctl_warn( 'The wave modified ocean stress components are used (ln_tauw=.true.) ', & 
     463                                  'This will override any other specification of the ocean stress' )  
     464         ENDIF 
     465      ELSE 
     466         IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_tauw .OR. ln_stcor .OR. ln_phioc .OR. ln_rough .OR. ln_zdfqiao ) &   
     467            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &  
     468            &                  'with drag coefficient (ln_cdgw =T) '  ,                        &  
     469            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &  
     470            &                  'or Stokes-Coriolis term (ln_stcor=T)',                         & 
     471            &                  'or ocean stress modification due to waves (ln_tauoc=T) ',      &  
     472            &                  'or ocean stress components from waves (ln_tauw=T) ',          &  
     473            &                  'or wave to ocean energy modification (ln_phioc=T) ',           &  
     474            &                  'or wave surface roughness (ln_rough=T) ',                      &  
     475            &                  'or Qiao vertical mixing formulation (ln_zdfqiao=T) ' ) 
     476      ENDIF  
     477      ! 
     478      IF( ln_cdgw ) THEN 
     479         IF( .NOT. cpl_wdrag ) THEN 
    96480            ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    97             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     481            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_cd structure' ) 
    98482            ! 
    99483                                   ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
    100484            IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
    101             CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
    102             ALLOCATE( cdn_wave(jpi,jpj) ) 
    103             cdn_wave(:,:) = 0.0 
    104         ENDIF 
    105          IF ( ln_sdw ) THEN 
    106             slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn 
    107             ALLOCATE( sf_sd(3), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    108             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     485            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 
     486         ENDIF 
     487         ALLOCATE( cdn_wave(jpi,jpj) ) 
     488      ENDIF 
     489 
     490      IF( ln_tauoc ) THEN 
     491         IF( .NOT. cpl_tauoc ) THEN 
     492            ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc 
     493            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauoc structure' ) 
    109494            ! 
    110             DO ifpr= 1, jpfld 
    111                ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
    112                IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
    113             END DO 
    114             CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
    115             ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj),uwavenum(jpi,jpj),vwavenum(jpi,jpj) ) 
    116             ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) ) 
    117             usd2d(:,:) = 0.0 ;  vsd2d(:,:) = 0.0 ; uwavenum(:,:) = 0.0 ; vwavenum(:,:) = 0.0 
    118             usd3d(:,:,:) = 0.0 ;vsd3d(:,:,:) = 0.0 ; wsd3d(:,:,:) = 0.0 
    119          ENDIF 
    120       ENDIF 
     495                                    ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   ) 
     496            IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 
     497            CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 
     498         ENDIF 
     499         ALLOCATE( tauoc_wave(jpi,jpj) ) 
     500      ENDIF 
     501 
     502      IF( ln_tauw ) THEN 
     503         IF( .NOT. cpl_tauw ) THEN 
     504            ALLOCATE( sf_tauw(2), STAT=ierror )           !* allocate and fill sf_wave with sn_tauwx/y 
     505            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauw structure' ) 
     506            ! 
     507            ALLOCATE( slf_j(2) ) 
     508            slf_j(1) = sn_tauwx 
     509            slf_j(2) = sn_tauwy 
     510                                    ALLOCATE( sf_tauw(1)%fnow(jpi,jpj,1)   ) 
     511                                    ALLOCATE( sf_tauw(2)%fnow(jpi,jpj,1)   ) 
     512            IF( slf_j(1)%ln_tint )  ALLOCATE( sf_tauw(1)%fdta(jpi,jpj,1,2) ) 
     513            IF( slf_j(2)%ln_tint )  ALLOCATE( sf_tauw(2)%fdta(jpi,jpj,1,2) ) 
     514            CALL fld_fill( sf_tauw, (/ slf_j /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 
     515         ENDIF 
     516         ALLOCATE( tauw_x(jpi,jpj) ) 
     517         ALLOCATE( tauw_y(jpi,jpj) ) 
     518      ENDIF 
     519 
     520      IF( ln_phioc ) THEN 
     521         IF( .NOT. cpl_phioc ) THEN 
     522            ALLOCATE( sf_phioc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_phioc 
     523            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_phioc structure' ) 
     524            ! 
     525                                    ALLOCATE( sf_phioc(1)%fnow(jpi,jpj,1)   ) 
     526            IF( sn_phioc%ln_tint )  ALLOCATE( sf_phioc(1)%fdta(jpi,jpj,1,2) ) 
     527            CALL fld_fill( sf_phioc, (/ sn_phioc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 
     528         ENDIF 
     529         ALLOCATE( rn_crban(jpi,jpj) ) 
     530      ENDIF 
     531 
     532      ! Find out how many fields have to be read from file if not coupled 
     533      jpfld=0 
     534      jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0   ;   jp_wfr=0 
     535      IF( ln_sdw ) THEN 
     536         IF( .NOT. cpl_sdrft ) THEN 
     537            jpfld  = jpfld + 1 
     538            jp_usd = jpfld 
     539            jpfld  = jpfld + 1 
     540            jp_vsd = jpfld 
     541         ENDIF 
     542         IF( .NOT. cpl_hsig .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 
     543            jpfld  = jpfld + 1 
     544            jp_hsw = jpfld 
     545         ENDIF 
     546         IF( .NOT. cpl_wper .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 
     547            jpfld  = jpfld + 1 
     548            jp_wmp = jpfld 
     549         ENDIF 
     550         IF( .NOT. cpl_wfreq .AND. nn_sdrift==jp_peakph ) THEN 
     551            jpfld  = jpfld + 1 
     552            jp_wfr = jpfld 
     553         ENDIF 
     554      ENDIF 
     555 
     556      IF( ln_rough .AND. .NOT. cpl_hsig .AND. jp_hsw==0 ) THEN 
     557         jpfld  = jpfld + 1 
     558         jp_hsw = jpfld 
     559      ENDIF 
     560 
     561      ! Read from file only the non-coupled fields  
     562      IF( jpfld > 0 ) THEN 
     563         ALLOCATE( slf_i(jpfld) ) 
     564         IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd 
     565         IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd 
     566         IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw 
     567         IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp 
     568         IF( jp_wfr > 0 )   slf_i(jp_wfr) = sn_wfr 
     569         ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift 
     570         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_sd structure' ) 
    121571         ! 
     572         DO ifpr= 1, jpfld 
     573            ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
     574            IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
     575         END DO 
    122576         ! 
    123       IF ( ln_cdgw ) THEN 
    124          CALL fld_read( kt, nn_fsbc, sf_cd )      !* read drag coefficient from external forcing 
    125          cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
    126       ENDIF 
    127       IF ( ln_sdw )  THEN 
    128           CALL fld_read( kt, nn_fsbc, sf_sd )      !* read drag coefficient from external forcing 
    129  
    130          ! Interpolate wavenumber, stokes drift into the grid_V and grid_V 
    131          !------------------------------------------------- 
    132  
    133          DO jj = 1, jpjm1 
    134             DO ji = 1, jpim1 
    135                uwavenum(ji,jj)=0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    136                &                                + sf_sd(3)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) ) 
    137  
    138                vwavenum(ji,jj)=0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    139                &                                + sf_sd(3)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) ) 
    140  
    141                usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(1)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    142                &                                + sf_sd(1)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) ) 
    143  
    144                vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(2)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    145                &                                + sf_sd(2)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) ) 
    146             END DO 
    147          END DO 
    148  
    149           !Computation of the 3d Stokes Drift 
    150           DO jk = 1, jpk 
    151              DO jj = 1, jpj-1 
    152                 DO ji = 1, jpi-1 
    153                    usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk)))) 
    154                    vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk)))) 
    155                 END DO 
    156              END DO 
    157              usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) ) 
    158              vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,jpj,jk)) ) 
    159           END DO 
    160  
    161           CALL wrk_alloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 
    162            
    163           udummy(:,:,:)=un(:,:,:) 
    164           vdummy(:,:,:)=vn(:,:,:) 
    165           hdivdummy(:,:,:)=hdivn(:,:,:) 
    166           rotdummy(:,:,:)=rotn(:,:,:) 
    167           un(:,:,:)=usd3d(:,:,:) 
    168           vn(:,:,:)=vsd3d(:,:,:) 
    169           CALL div_cur(kt) 
    170       !                                           !------------------------------! 
    171       !                                           !     Now Vertical Velocity    ! 
    172       !                                           !------------------------------! 
    173           z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
    174  
    175           z1_2dt = 1.e0 / z2dt 
    176           DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence 
    177              ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 
    178              wsd3d(:,:,jk) = wsd3d(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        & 
    179                 &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    & 
    180                 &                         * tmask(:,:,jk) * z1_2dt 
    181 #if defined key_bdy 
    182              wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 
    183 #endif 
    184           END DO 
    185           hdivn(:,:,:)=hdivdummy(:,:,:) 
    186           rotn(:,:,:)=rotdummy(:,:,:) 
    187           vn(:,:,:)=vdummy(:,:,:) 
    188           un(:,:,:)=udummy(:,:,:) 
    189           CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 
    190       ENDIF 
    191    END SUBROUTINE sbc_wave 
    192        
     577         CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 
     578      ENDIF 
     579 
     580      IF( ln_sdw ) THEN 
     581         ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) ) 
     582         ALLOCATE( wmp  (jpi,jpj)  ) 
     583         ALLOCATE( wfreq (jpi,jpj) ) 
     584         ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     ) 
     585         ALLOCATE( div_sd(jpi,jpj) ) 
     586         ALLOCATE( tsd2d (jpi,jpj) ) 
     587         usd(:,:,:) = 0._wp 
     588         vsd(:,:,:) = 0._wp 
     589         wsd(:,:,:) = 0._wp 
     590         ! Wave number needed only if ln_zdfqiao=T 
     591         IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN 
     592            ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum 
     593            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wn structure' ) 
     594                                   ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
     595            IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
     596            CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'read wave input', 'namsbc_wave' ) 
     597         ENDIF 
     598         ALLOCATE( wnum(jpi,jpj) ) 
     599      ENDIF 
     600 
     601      IF( ln_sdw .OR. ln_rough ) THEN 
     602         ALLOCATE( hsw  (jpi,jpj) ) 
     603      ENDIF 
     604      ! 
     605   END SUBROUTINE sbc_wave_init 
     606 
    193607   !!====================================================================== 
    194608END MODULE sbcwave 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r8058 r11286  
    66   !! History :  2.0  !  2005-11  (G. Madec)  Original code 
    77   !!            3.3  !  2010-09  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport 
     8   !!            3.6  !  2015-06  (E. Clementi) Addition of Stokes drift in case of wave coupling 
    89   !!            4.0  !  2011-06  (G. Madec)  Addition of Mixed Layer Eddy parameterisation 
    910   !!---------------------------------------------------------------------- 
    … …  
    3435   USE timing          ! Timing 
    3536   USE sbc_oce 
    36    USE diaptr          ! Poleward heat transport  
    37  
     37   USE sbcwave        ! wave module   
     38   USE sbc_oce        ! surface boundary condition: ocean   
     39   USE diaptr         ! Poleward heat transport 
    3840 
    3941   IMPLICIT NONE 
    … …  
    8587      CALL wrk_alloc( jpi, jpj, jpk, zun, zvn, zwn ) 
    8688      !                                          ! set time step 
     89      zun(:,:,:) = 0.0   
     90      zvn(:,:,:) = 0.0   
     91      zwn(:,:,:) = 0.0   
     92      ! 
    8793      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000 
    8894         r2dtra(:) =  rdttra(:)                          ! = rdtra (restarting with Euler time stepping) 
    … …  
    94100      ! 
    95101      !                                               !==  effective transport  ==! 
    96       DO jk = 1, jpkm1 
    97          zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only 
    98          zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
    99          zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      * wn(:,:,jk) 
    100       END DO 
     102      IF(ln_wave .AND. ln_sdw) THEN   
     103         DO jk = 1, jpkm1                                                                      ! eulerian transport + Stokes Drift   
     104            zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * ( un(:,:,jk) + usd(:,:,jk) )  
     105            zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * ( vn(:,:,jk) + vsd(:,:,jk) )   
     106            zwn(:,:,jk) = e1e2t(:,:)               * ( wn(:,:,jk) + wsd(:,:,jk) )   
     107         END DO   
     108      ELSE   
     109         DO jk = 1, jpkm1   
     110            zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)             ! eulerian transport only   
     111            zvn(:,:,jk) = e1v  (:,:) * fse3v(:,:,jk) * vn(:,:,jk)   
     112            zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk)   
     113         END DO   
     114      ENDIF 
    101115      ! 
    102116      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90

    r8058 r11286  
    1010   USE in_out_manager ! I/O manager 
    1111   USE lib_mpp        ! MPP library 
     12   USE sbc_oce, ONLY : ln_zdfqiao 
    1213 
    1314   IMPLICIT NONE 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r8058 r11286  
    2424   USE phycst         ! physical constants 
    2525   USE zdfmxl         ! mixed layer 
     26   USE sbcwave 
     27   ! 
    2628   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    2729   USE lib_mpp        ! MPP manager 
    … …  
    4648   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points 
    4749   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustarb2 !: Squared bottom  velocity scale at T-points 
     50 
     51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rsbc_tke1  
     52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rsbc_tke3  
     53   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rsbc_psi1  
     54   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rtrans 
    4855 
    4956   !                              !! ** Namelist  namzdf_gls  ** 
    … …  
    5966   REAL(wp) ::   rn_emin           ! minimum value of TKE (m2/s2) 
    6067   REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 
    61    REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing 
     68   REAL(wp) ::   rn_crban_default  ! Craig and Banner constant for surface breaking waves mixing 
    6269   REAL(wp) ::   rn_hsro           ! Minimum surface roughness 
    63    REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1)  
     70   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met = 1, 3)  
    6471 
    6572   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters 
    … …  
    9097   REAL(wp) ::   rm7           =  0.0_wp 
    9198   REAL(wp) ::   rm8           =  0.318_wp 
    92    REAL(wp) ::   rtrans        =  0.1_wp 
     99   REAL(wp) ::   rtrans_default =  0.1_wp 
    93100   REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters 
    94    REAL(wp) ::   rsbc_tke1, rsbc_tke2, rfact_tke                  !     -           -           -        - 
    95    REAL(wp) ::   rsbc_psi1, rsbc_psi2, rfact_psi                  !     -           -           -        - 
     101   REAL(wp) ::   rsbc_tke1_default, rsbc_tke2, rfact_tke          !     -           -           -        -  
     102   REAL(wp) ::   rsbc_psi2, rsbc_psi3, rfact_psi                  !     -           -           -        - 
    96103   REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        - 
    97104   REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        - 
    … …  
    116123      !!---------------------------------------------------------------------- 
    117124      ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     & 
    118          &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc ) 
     125         &      rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) ,     &  
     126         &      rsbc_psi1(jpi,jpj), rtrans(jpi,jpj),         & 
     127         &      ustars2(jpi,jpj), ustarb2(jpi,jpj)                      , STAT= zdf_gls_alloc ) 
    119128         ! 
    120129      IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc ) 
    … …  
    161170      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp 
    162171 
     172      ! variable initialization  
     173      IF( ln_phioc ) THEN  
     174         IF( nn_wmix==jp_janssen ) THEN  
     175            rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.Dirichlet + Wave breaking   
     176            rsbc_tke3(:,:) = rdt * rn_crban(:,:)                                           ! Neumann + Wave breaking  
     177            rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn                   ! Dirichlet + Wave breaking  
     178         ELSE IF( nn_wmix==jp_craigbanner ) THEN  
     179            rsbc_tke1(:,:) = -3._wp/2._wp*rn_crban(:,:)*ra_sf*rl_sf  
     180            rsbc_tke3(:,:) = rdt * rn_crban(:,:) / rl_sf  
     181            rtrans(:,:) = 0.2_wp / MAX( rsmall, rsbc_tke1(:,:)**(1./(-ra_sf*3._wp/2._wp))-1._wp )  
     182         ENDIF  
     183      ENDIF  
     184 
    163185      IF( kt /= nit000 ) THEN   ! restore before value to compute tke 
    164186         avt (:,:,:) = avt_k (:,:,:) 
    … …  
    197219         zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10) 
    198220         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 
    199       ! 
     221      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)   
     222         zhsro = MAX(rn_frac_hs * hsw, rn_hsro) 
    200223      END SELECT 
    201224 
    … …  
    311334      ! --------------------------------- 
    312335      ! 
     336      IF( ln_phioc ) THEN  
     337         SELECT CASE ( nn_bc_surf )  
     338         !  
     339         CASE ( 0 )             ! Dirichlet case  
     340            IF( nn_wmix==jp_janssen ) THEN  
     341               ! First level  
     342               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin )  
     343               z_elem_a(:,:,1) = en(:,:,1)  
     344               z_elem_c(:,:,1) = 0._wp  
     345               z_elem_b(:,:,1) = 1._wp  
     346  
     347               ! One level below  
     348               en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin )  
     349               z_elem_a(:,:,2) = 0._wp  
     350               z_elem_c(:,:,2) = 0._wp  
     351               z_elem_b(:,:,2) = 1._wp  
     352            ELSE IF( nn_wmix==jp_craigbanner ) THEN  
     353               en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp)  
     354               en(:,:,1) = MAX(en(:,:,1), rn_emin)   
     355               z_elem_a(:,:,1) = en(:,:,1)  
     356               z_elem_c(:,:,1) = 0._wp  
     357               z_elem_b(:,:,1) = 1._wp  
     358               !   
     359               ! One level below  
     360               en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:) * ((zhsro(:,:)+fsdepw(:,:,2)) &  
     361                   &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)  
     362               en(:,:,2) = MAX(en(:,:,2), rn_emin )  
     363               z_elem_a(:,:,2) = 0._wp   
     364               z_elem_c(:,:,2) = 0._wp  
     365               z_elem_b(:,:,2) = 1._wp  
     366               !  
     367            ENDIF  
     368         CASE ( 1 )             ! Neumann boundary condition on d(e)/dz  
     369            IF( nn_wmix==jp_janssen ) THEN  
     370               ! Dirichlet conditions at k=1  
     371               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin )  
     372               z_elem_a(:,:,1) = en(:,:,1)  
     373               z_elem_c(:,:,1) = 0._wp  
     374               z_elem_b(:,:,1) = 1._wp  
     375               ! at k=2, set de/dz=Fw  
     376               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b  
     377               z_elem_a(:,:,2) = 0._wp  
     378               zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf)  
     379               en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2)  
     380            ELSE IF( nn_wmix==jp_craigbanner ) THEN  
     381               ! Dirichlet conditions at k=1  
     382               en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp)  
     383               en(:,:,1)       = MAX(en(:,:,1), rn_emin)        
     384               z_elem_a(:,:,1) = en(:,:,1)  
     385               z_elem_c(:,:,1) = 0._wp  
     386               z_elem_b(:,:,1) = 1._wp  
     387               !  
     388               ! at k=2, set de/dz=Fw  
     389               !cbr  
     390               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b  
     391               z_elem_a(:,:,2) = 0._wp  
     392               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:)) ))  
     393               zflxs(:,:)      = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * zkar(:,:) &  
     394                    &                           * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)  
     395 
     396               en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2)  
     397               !  
     398            ENDIF  
     399         END SELECT  
     400      ELSE 
    313401      SELECT CASE ( nn_bc_surf ) 
    314402      ! 
    315403      CASE ( 0 )             ! Dirichlet case 
    316404      ! First level 
    317       en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
     405      en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 
    318406      en(:,:,1) = MAX(en(:,:,1), rn_emin)  
    319407      z_elem_a(:,:,1) = en(:,:,1) 
    … …  
    322410      !  
    323411      ! One level below 
    324       en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 
    325           &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
     412      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
    326413      en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
    327414      z_elem_a(:,:,2) = 0._wp  
    … …  
    331418      ! 
    332419      CASE ( 1 )             ! Neumann boundary condition on d(e)/dz 
    333       ! 
    334420      ! Dirichlet conditions at k=1 
    335       en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
     421      en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 
    336422      en(:,:,1)       = MAX(en(:,:,1), rn_emin)       
    337423      z_elem_a(:,:,1) = en(:,:,1) 
    … …  
    343429      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    344430      z_elem_a(:,:,2) = 0._wp 
    345       zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 
    346       zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 
    347            &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
     431      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:)) )) 
     432      zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
    348433 
    349434      en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 
    … …  
    351436      ! 
    352437      END SELECT 
     438      ENDIF ! ln_phioc 
    353439 
    354440      ! Bottom boundary condition on tke 
    … …  
    536622      ! 
    537623      CASE ( 0 )             ! Dirichlet boundary conditions 
    538       ! 
    539       ! Surface value 
    540       zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
    541       psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    542       z_elem_a(:,:,1) = psi(:,:,1) 
    543       z_elem_c(:,:,1) = 0._wp 
    544       z_elem_b(:,:,1) = 1._wp 
    545       ! 
    546       ! One level below 
    547       zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 
    548       zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 
    549       psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    550       z_elem_a(:,:,2) = 0._wp 
    551       z_elem_c(:,:,2) = 0._wp 
    552       z_elem_b(:,:,2) = 1._wp 
    553       !  
    554       ! 
     624         IF( ln_phioc ) THEN   ! Wave induced mixing case  
     625                               ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2  
     626                               ! balance between the production and the  
     627                               ! dissipation terms including the wave effect  
     628            IF( nn_wmix==jp_janssen ) THEN  
     629               ! Surface value  
     630               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic  
     631               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     632               z_elem_a(:,:,1) = psi(:,:,1)        
     633               z_elem_c(:,:,1) = 0._wp  
     634               z_elem_b(:,:,1) = 1._wp  
     635               !  
     636               ! One level below  
     637               zex1 = (rmm*ra_sf+rnn)  
     638               zex2 = (rmm*ra_sf)  
     639               zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2  
     640               psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1)  
     641               z_elem_a(:,:,2) = 0._wp  
     642               z_elem_c(:,:,2) = 0._wp  
     643               z_elem_b(:,:,2) = 1._wp  
     644               !   
     645               !  
     646            ELSE IF( nn_wmix==jp_craigbanner ) THEN  
     647               ! Surface value  
     648               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic  
     649               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     650               z_elem_a(:,:,1) = psi(:,:,1)  
     651               z_elem_c(:,:,1) = 0._wp  
     652               z_elem_b(:,:,1) = 1._wp  
     653               !  
     654               ! One level below  
     655               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdepw(:,:,2)/zhsro(:,:) )))  
     656               zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)  
     657               psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     658               z_elem_a(:,:,2) = 0._wp  
     659               z_elem_c(:,:,2) = 0._wp  
     660               z_elem_b(:,:,2) = 1._wp  
     661               !   
     662               !  
     663            ENDIF  
     664         ELSE                  ! Wave induced mixing case with default values  
     665            ! Surface value  
     666            zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic  
     667            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     668            z_elem_a(:,:,1) = psi(:,:,1)  
     669            z_elem_c(:,:,1) = 0._wp  
     670            z_elem_b(:,:,1) = 1._wp  
     671            !  
     672            ! One level below  
     673            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdepw(:,:,2)/zhsro(:,:) )))  
     674            zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)  
     675            psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     676            z_elem_a(:,:,2) = 0._wp  
     677            z_elem_c(:,:,2) = 0._wp  
     678            z_elem_b(:,:,2) = 1._wp  
     679            !   
     680            !  
     681         ENDIF 
    555682      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz 
    556       ! 
    557       ! Surface value: Dirichlet 
    558       zdep(:,:)       = zhsro(:,:) * rl_sf 
    559       psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    560       z_elem_a(:,:,1) = psi(:,:,1) 
    561       z_elem_c(:,:,1) = 0._wp 
    562       z_elem_b(:,:,1) = 1._wp 
    563       ! 
    564       ! Neumann condition at k=2 
    565       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    566       z_elem_a(:,:,2) = 0._wp 
    567       ! 
    568       ! Set psi vertical flux at the surface: 
    569       zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
    570       zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
    571       zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
    572       zdep(:,:) =  rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 
    573              & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 
    574       zflxs(:,:) = zdep(:,:) * zflxs(:,:) 
    575       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
    576  
    577       !    
    578       ! 
     683         IF( ln_phioc ) THEN  ! Wave induced mixing case with forced/coupled fields  
     684            IF( nn_wmix==jp_janssen ) THEN  
     685               !  
     686               zdep(:,:) = rl_sf * zhsro(:,:)  
     687               psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     688               z_elem_a(:,:,1) = psi(:,:,1)  
     689               z_elem_c(:,:,1) = 0._wp  
     690               z_elem_b(:,:,1) = 1._wp  
     691               !  
     692               ! Neumann condition at k=2  
     693               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b  
     694               z_elem_a(:,:,2) = 0._wp  
     695               !  
     696               ! Set psi vertical flux at the surface:  
     697               zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf)  
     698               zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) &  
     699                  &                   * en(:,:,1)**rmm * zdep  
     700               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)  
     701               !  
     702            ELSE IF( nn_wmix==jp_craigbanner ) THEN  
     703               ! Surface value: Dirichlet  
     704               zdep(:,:)       = zhsro(:,:) * rl_sf  
     705               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     706               z_elem_a(:,:,1) = psi(:,:,1)  
     707               z_elem_c(:,:,1) = 0._wp  
     708               z_elem_b(:,:,1) = 1._wp  
     709               !  
     710               ! Neumann condition at k=2  
     711               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b  
     712               z_elem_a(:,:,2) = 0._wp  
     713               !  
     714               ! Set psi vertical flux at the surface:  
     715               zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope  
     716               zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)  
     717               zflxs(:,:) = (rnn + rsbc_tke1(:,:) * (rnn + rmm*ra_sf) * zdep(:,:)) * &  
     718                          (1._wp + rsbc_tke1(:,:) * zdep(:,:))**(2._wp*rmm/3._wp-1_wp)  
     719               zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &  
     720                      & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.)  
     721               zflxs(:,:) = zdep(:,:) * zflxs(:,:)  
     722               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)  
     723               !     
     724            ENDIF  
     725         ELSE                 ! Wave induced mixing case with default values  
     726            ! Surface value: Dirichlet  
     727            zdep(:,:)       = zhsro(:,:) * rl_sf  
     728            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     729            z_elem_a(:,:,1) = psi(:,:,1)  
     730            z_elem_c(:,:,1) = 0._wp  
     731            z_elem_b(:,:,1) = 1._wp  
     732            !  
     733            ! Neumann condition at k=2  
     734            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b  
     735            z_elem_a(:,:,2) = 0._wp  
     736            !  
     737            ! Set psi vertical flux at the surface:  
     738            zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope  
     739            zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)  
     740            zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * &  
     741                       (1._wp + rsbc_tke1_default * zdep(:,:))**(2._wp*rmm/3._wp-1_wp)  
     742            zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &  
     743                   & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.)  
     744            zflxs(:,:) = zdep(:,:) * zflxs(:,:)  
     745            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)  
     746            !     
     747            !  
     748         ENDIF 
    579749      END SELECT 
    580750 
    … …  
    8671037      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 
    8681038         &            rn_clim_galp, ln_sigpsi, rn_hsro,      & 
    869          &            rn_crban, rn_charn, rn_frac_hs,        & 
     1039         &            rn_crban_default, rn_charn, rn_frac_hs,& 
    8701040         &            nn_bc_surf, nn_bc_bot, nn_z0_met,      & 
    8711041         &            nn_stab_func, nn_clos 
    … …  
    8951065         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot 
    8961066         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi 
    897          WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban 
     1067         WRITE(numout,*) '      Craig and Banner coefficient (default)        rn_crban       = ', rn_crban_default 
    8981068         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn 
    8991069         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met 
    … …  
    9091079 
    9101080      !                                !* Check of some namelist values 
    911       IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 
    912       IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 
    913       IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' ) 
    914       IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 
    915       IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) 
     1081      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )   
     1082      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )   
     1083      IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' )   
     1084      IF( nn_z0_met == 3 .AND. .NOT.ln_rough ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_rough=T' )   
     1085      IF( nn_z0_met .NE. 3 .AND. ln_rough ) THEN  
     1086         CALL ctl_warn('W A R N I N G:  ln_rough=.TRUE. but nn_z0_met is not 3 - resetting nn_z0_met to 3')   
     1087         nn_z0_met = 3   
     1088      ENDIF 
     1089      IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' )   
     1090      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 
    9161091 
    9171092      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure 
    … …  
    10851260               &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 
    10861261 
    1087       IF ( rn_crban==0._wp ) THEN 
     1262      IF( .NOT. ln_phioc .AND. rn_crban_default==0._wp ) THEN 
    10881263         rl_sf = vkarmn 
    10891264      ELSE 
    … …  
    11241299      rc03  = rc02 * rc0 
    11251300      rc04  = rc03 * rc0 
    1126       rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking 
    1127       rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking  
    1128       zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
    1129       rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer  
     1301      rsbc_tke1_default = -3._wp/2._wp*rn_crban_default*ra_sf*rl_sf  
     1302      rsbc_tke2         = rdt * rn_crban_default / rl_sf  
     1303      zcr               = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp )  
     1304      rtrans_default    = 0.2_wp / zcr 
    11301305      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness 
    11311306      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness  

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90

    r8058 r11286  
    5353      INTEGER ::   ios 
    5454      !! 
    55       NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,   & 
    56          &              ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp 
     55      NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,  &   
     56         &        ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp 
    5757      !!---------------------------------------------------------------------- 
    5858 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r10728 r11286  
    6969   USE icbini          ! handle bergs, initialisation 
    7070   USE icbstp          ! handle bergs, calving, themodynamics and transport 
     71   USE sbccpl 
    7172   USE cpl_oasis3      ! OASIS3 coupling 
    7273   USE c1d             ! 1D configuration 
    … …  
    175176            CALL stp                         ! AGRIF: time stepping 
    176177#else 
    177             CALL stp( istp )                 ! standard time stepping 
     178            IF (lk_oasis) CALL sbc_cpl_snd( istp )  ! Coupling to atmos    
     179            CALL stp( istp )    
     180            ! We don't couple on the final timestep because    
     181            ! our restart file has already been written    
     182            ! and contains all the necessary data for a    
     183            ! restart. sbc_cpl_snd could be called here    
     184            ! but it would require    
     185            ! a) A test to ensure it was not performed    
     186            !    on the very last time-step    
     187            ! b) the presence of another call to    
     188            !    sbc_cpl_snd call prior to the main DO loop    
     189            ! This solution produces identical results    
     190            ! with fewer lines of code.   
    178191#endif 
    179192            istp = istp + 1 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/step.F90

    r10728 r11286  
    2525   !!            3.4  !  2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal 
    2626   !!                 !  2012-07  (J. Simeon, G. Madec, C. Ethe) Online coarsening of outputs 
     27   !!            3.6  !  2014-10  (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves 
    2728   !!            3.7  !  2014-04  (F. Roquet, G. Madec) New equations of state 
    2829   !!---------------------------------------------------------------------- 
    … …  
    7273      !!              -8- Outputs and diagnostics 
    7374      !!---------------------------------------------------------------------- 
    74       INTEGER ::   jk      ! dummy loop indice 
     75      INTEGER ::   ji,jj,jk ! dummy loop indice 
    7576      INTEGER ::   indic    ! error indicator if < 0 
    7677      INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt) 
    … …  
    105106      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    106107      IF( lk_tide    )   CALL sbc_tide( kstp ) 
    107       IF( lk_bdy     )  THEN 
    108          IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib  
    109                          CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries 
    110       ENDIF 
    111                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice) 
     108                         CALL sbc     ( kstp )        ! Sea Boundary Condition (including sea-ice) 
    112109                                                      ! clem: moved here for bdy ice purpose 
    113110      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    132129      IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz 
    133130      IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz 
     131      IF( ln_zdfqiao )   CALL zdf_qiao( kstp )           ! Qiao vertical mixing    
     132      ! 
    134133      IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz 
    135134      IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value) 
    … …  
    399398                               CALL ctl_stop( 'step: indic < 0' ) 
    400399                               CALL dia_wri_state( 'output.abort', kstp ) 
     400                               CALL ctl_stop('STOP','NEMO failure in stp') 
    401401      ENDIF 
    402402      IF( ln_harm_ana_store   )   CALL harm_ana( kstp )        ! Harmonic analysis of tides  
    … …  
    411411      ! Coupled mode 
    412412      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    413       IF( lk_oasis         )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges 
     413      !IF( lk_oasis         )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges 
    414414      ! 
    415415#if defined key_iomput 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r8561 r11286  
    2929   USE sbctide          ! Tide initialisation 
    3030   USE sbcapr           ! surface boundary condition: ssh_ib required by bdydta  
     31   USE sbcwave          ! Wave intialisation 
    3132 
    3233   USE traqsr           ! solar radiation penetration      (tra_qsr routine) 
    … …  
    8687   USE zdfric           ! Richardson vertical mixing       (zdf_ric routine) 
    8788   USE zdfmxl           ! Mixed-layer depth                (zdf_mxl routine) 
     89   USE zdfqiao          !Qiao module wave induced mixing   (zdf_qiao routine) 
    8890 
    8991   USE zpshde           ! partial step: hor. derivative     (zps_hde routine) 

  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_xeps/NEMOGCM/SETTE/sette.sh

    r5588 r11286  
    8888# 
    8989# Compiler among those in NEMOGCM/ARCH 
    90 COMPILER=X64_ADA 
    91 export BATCH_COMMAND_PAR="llsubmit" 
     90module load cray-netcdf-hdf5parallel   
     91COMPILER=XC40_METO   
     92export BATCH_COMMAND_PAR="qsub 
    9293export BATCH_COMMAND_SEQ=$BATCH_COMMAND_PAR 
    9394export INTERACT_FLAG="no" 
Note: See TracChangeset for help on using the changeset viewer.