Changeset 10478


Ignore:
Timestamp:
2019-01-09T12:11:55+01:00 (15 months ago)
Author:
jcastill
Message:

Merged branch UKMO/r6232_HZG_WAVE-coupling@9868

Location:
branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM
Files:
16 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/CONFIG/SHARED/namelist_ref

    r10473 r10478  
    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 
     
    265266                           !     =2 annual global mean of e-p-r set to zero 
    266267   ln_wave = .false.       !  Activate coupling with wave (T => fill namsbc_wave)    
    267    ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => ln_wave=.true. & fill namsbc_wave)    
    268    ln_sdw  = .false.       !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave)     
    269    ln_tauoc= .false.       !  Activate ocean stress modified by external wave induced stress (T => ln_wave=.true. & fill namsbc_wave)    
    270    ln_stcor= .false.       !  Activate Stokes Coriolis term (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave 
    271268   nn_lsm  = 0             !  =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) , 
    272269                           !  =1:n number of iterations of land/sea mask application for input fields (fill land/sea mask filename field) 
     
    276273                           !  = 1  Average and redistribute per-category fluxes, forced mode only, not yet implemented coupled 
    277274                           !  = 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 
    278280/ 
    279281!----------------------------------------------------------------------- 
     
    298300   sn_emp      = 'emp'       ,        24         , 'emp'     , .false.      , .false., 'yearly'  , ''       , ''       , '' 
    299301 
    300    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 
    301306/ 
    302307!----------------------------------------------------------------------- 
     
    383388   sn_rcv_mslp   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
    384389   sn_rcv_phioc  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
    385    sn_rcv_sdrfx  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
    386    sn_rcv_sdrfy  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
     390   sn_rcv_sdrft  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
    387391   sn_rcv_wper   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
     392   sn_rcv_wfreq  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
    388393   sn_rcv_wnum   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
    389    sn_rcv_wstrf  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
     394   sn_rcv_tauoc  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''   
     395   sn_rcv_tauw   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
    390396   sn_rcv_wdrag  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
    391397! 
     
    936942   ln_zdfexp   = .false.   !  time-stepping: split-explicit (T) or implicit (F) time stepping 
    937943   nn_zdfexp   =    3            !  number of sub-timestep for ln_zdfexp=T 
    938    ln_zdfqiao  = .false.   !  Enhanced wave vertical mixing Qiao (2010) (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave) 
    939944/ 
    940945!----------------------------------------------------------------------- 
     
    10001005   rn_clim_galp  = 0.267   !  galperin limit 
    10011006   ln_sigpsi     = .true.  !  Activate or not Burchard 2001 mods on psi schmidt number in the wb case 
    1002    rn_crban      = 100.    !  Craig and Banner 1994 constant for wb tke flux 
     1007   rn_crban_default = 100.    !  Craig and Banner 1994 constant for wb tke flux 
    10031008   rn_charn      = 70000.  !  Charnock constant for wb induced roughness length 
    10041009   rn_hsro       =  0.02   !  Minimum surface roughness 
     
    12941299!              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
    12951300!              !             !  (if <0  months)  !   name       !   (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      ! 
     1301   ln_cdgw     = .false.   !  Neutral drag coefficient read from wave model  
     1302   ln_sdw      = .false.   !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift  
     1303   ln_stcor    = .false.   !  Activate Stokes Coriolis term  
     1304   ln_tauoc    = .false.   !  Activate ocean stress modified by external wave induced stress  
     1305   ln_tauw     = .false.   !  Activate ocean stress components from wave model  
     1306   ln_phioc    = .false.   !  Activate wave to ocean energy  
     1307   ln_rough    = .false.   !  Wave roughness equals the significant wave height  
     1308   ln_zdfqiao  = .false.   !  Enhanced wave vertical mixing Qiao (2010)  
     1309   nn_sdrift   =  0   !  Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift  
     1310                           ! = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]  
     1311                           ! = 1 Phillips:                      v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))]  
     1312   nn_wmix     =  0   !  type of wave breaking mixing  
     1313                           ! = 0 Craig and Banner formulation (original NEMO formulation)   
     1314                           ! = 1 Janssen formulation (no assumption of direct energy conversion) 
    12961315   sn_cdg      =  'cdw_wave' ,        1          , 'drag_coeff' ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    12971316   sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
     
    12991318   sn_hsw      =  'sdw_wave' ,        1          , 'hs'         ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
    13001319   sn_wmp      =  'sdw_wave' ,        1          , 'wmp'        ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
     1320   sn_wfr      =  'sdw_wave' ,        1          , 'wave_freq'  ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    13011321   sn_wnum     =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
    13021322   sn_tauoc    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
     1323   sn_tauwx    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
     1324   sn_tauwy    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
     1325   sn_phioc    =  'sdw_wave' ,        1          , 'wave_energy',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''   
     1326   cn_dir      = './'  !  root directory for the location of wave forcing files 
    13031327 
    1304    cn_dir  = './'  !  root directory for the location of drag coefficient files 
    13051328/ 
    13061329!----------------------------------------------------------------------- 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r10473 r10478  
    10941094 
    10951095      ! 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 
    10961101      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature 
    10971102      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity 
     
    11111116      END IF 
    11121117        
    1113       IF( ln_wave .AND. ln_sdw ) THEN   
    1114          CALL histwrite( id_i, "sdzocrtx", kt, usd           , jpi*jpj*jpk, idex)     ! now StokesDrift i-velocity   
    1115          CALL histwrite( id_i, "sdmecrty", kt, vsd           , jpi*jpj*jpk, idex)     ! now StokesDrift j-velocity   
    1116          CALL histwrite( id_i, "sdvecrtz", kt, wsd           , jpi*jpj*jpk, idex)     ! now StokesDrift k-velocity   
    1117       ENDIF   
    1118  
    11191118      ! 3. Close the file 
    11201119      ! ----------------- 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r10473 r10478  
    470470                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     ) 
    471471      ENDIF 
    472       !   
    473       IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary   
    474          zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)   
    475       ENDIF   
    476       !   
    477472#if defined key_asminc 
    478473      !                                         ! Include the IAU weighted SSH increment 
     
    488483      ! 
    489484      ! ----------------------------------------------------------------------- 
    490       !  Phase 2 : Integration of the barotropic equations  
     485      !  Phase 2 : Integration of the barotropic equations 
    491486      ! ----------------------------------------------------------------------- 
    492487      ! 
    493488      !                                             ! ==================== ! 
    494489      !                                             !    Initialisations   ! 
    495       !                                             ! ==================== !   
     490      !                                             ! ==================== ! 
    496491      ! Initialize barotropic variables:       
    497492      IF( ll_init )THEN 
     
    502497         ub_e   (:,:) = 0._wp 
    503498         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(:,:)   
    504504      ENDIF 
    505505      ! 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r10473 r10478  
    182182         ELSE 
    183183                             CALL vor_een( kt, ntot, un, vn, ua, va )    ! total vorticity  
    184             IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend 
     184            IF( ln_stcor )   CALL vor_een( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend 
    185185         ENDIF 
    186186         ! 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r10473 r10478  
    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 
     
    6466   LOGICAL , PUBLIC ::   ln_wave        !: true if some coupling with wave model 
    6567   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 
     68   LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d Stokes drift from wave model  
     69   LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used 
    6770   LOGICAL , PUBLIC ::   ln_tauoc       !: true if normalized stress from wave is used   
    68    LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used 
     71   LOGICAL , PUBLIC ::   ln_tauw        !: true if ocean stress components from wave is used   
     72   LOGICAL , PUBLIC ::   ln_phioc       !: true if wave energy to ocean is used   
     73   LOGICAL , PUBLIC ::   ln_rough       !: true if wave roughness equals significant wave height  
     74   LOGICAL , PUBLIC ::   ln_zdfqiao     !: Enhanced wave vertical mixing Qiao(2010) formulation flag  
     75   INTEGER , PUBLIC ::   nn_drag        ! type of formula to calculate wind stress from wind components  
     76   INTEGER , PUBLIC ::   nn_wmix        ! type of wave breaking mixing 
    6977   ! 
    7078   LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs 
     
    93101   INTEGER , PUBLIC, PARAMETER ::   jp_iam_sas  = 2      !: Multi executable configuration - SAS component 
    94102                                                         !  (internal OASIS coupling) 
     103   !!----------------------------------------------------------------------  
     104   !!           wind stress definition  
     105   !!----------------------------------------------------------------------  
     106   INTEGER, PUBLIC, PARAMETER ::   jp_ukmo  = 0        ! UKMO SHELF formulation  
     107   INTEGER, PUBLIC, PARAMETER ::   jp_std   = 1        ! standard formulation with forced or coupled drag coefficient   
     108   INTEGER, PUBLIC, PARAMETER ::   jp_const = 2        ! standard formulation with constant drag coefficient   
     109   INTEGER, PUBLIC, PARAMETER ::   jp_mcore = 3        ! momentum calculated from core forcing fields   
     110 
     111   !!----------------------------------------------------------------------  
     112   !!           Wave mixing vertical parameterization  
     113   !!----------------------------------------------------------------------  
     114   INTEGER, PUBLIC, PARAMETER ::   jp_craigbanner = 0  ! Craig and Banner formulation (original NEMO formulation -  
     115                                                       ! direct conversion of mechanical to turbulent energy)  
     116   INTEGER, PUBLIC, PARAMETER ::   jp_janssen     = 1  ! Janssen formulation - no assumption on direct energy conversion 
    95117   !!---------------------------------------------------------------------- 
    96118   !!              Ocean Surface Boundary Condition fields 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90

    r8058 r10478  
    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_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r10473 r10478  
    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_wave .AND. 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_wave .AND. 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_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r10473 r10478  
    116116   INTEGER, PARAMETER ::   jpr_wper   = 48            ! Mean wave period   
    117117   INTEGER, PARAMETER ::   jpr_wnum   = 49            ! Mean wavenumber   
    118    INTEGER, PARAMETER ::   jpr_wstrf  = 50            ! Stress fraction adsorbed by waves   
     118   INTEGER, PARAMETER ::   jpr_tauoc  = 50            ! Stress fraction adsorbed by waves 
    119119   INTEGER, PARAMETER ::   jpr_wdrag  = 51            ! Neutral surface drag coefficient   
    120    INTEGER, PARAMETER ::   jprcv      = 51            ! total number of fields received 
     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 
    121124 
    122125   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere 
     
    170173   TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev    
    171174   ! Received from waves    
    172    TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrfx,sn_rcv_sdrfy,sn_rcv_wper,sn_rcv_wnum,sn_rcv_wstrf,sn_rcv_wdrag 
     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 
    173178   ! Other namelist parameters                        ! 
    174179   INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 
     
    183188   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky) 
    184189     
    185    REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2]    
    186    REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0 
    187  
    188190   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument 
    189191 
     
    247249         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,      &    
    248250         &                  sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev  , sn_rcv_hsig  , sn_rcv_phioc ,   &    
    249          &                  sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_wstrf ,   &   
    250          &                  sn_rcv_wdrag, sn_rcv_qns  , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   ,   &   
    251          &                  sn_rcv_iceflx,sn_rcv_co2  , nn_cplmodel  , ln_usecplmask, sn_rcv_mslp 
     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 
    252255      !!--------------------------------------------------------------------- 
    253256      ! 
     
    290293         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 
    291294         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  ), ')' 
    292296         WRITE(numout,*)'      significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')'    
    293297         WRITE(numout,*)'      wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')'    
    294          WRITE(numout,*)'      Surface Stokes drift grid u     = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')'    
    295          WRITE(numout,*)'      Surface Stokes drift grid v     = ', TRIM(sn_rcv_sdrfy%cldes ), ' (', TRIM(sn_rcv_sdrfy%clcat ), ')'    
     298         WRITE(numout,*)'      Surface Stokes drift u,v        = ', TRIM(sn_rcv_sdrft%cldes ), ' (', TRIM(sn_rcv_sdrft%clcat ), ')' 
    296299         WRITE(numout,*)'      Mean wave period                = ', TRIM(sn_rcv_wper%cldes  ), ' (', TRIM(sn_rcv_wper%clcat  ), ')'    
    297300         WRITE(numout,*)'      Mean wave number                = ', TRIM(sn_rcv_wnum%cldes  ), ' (', TRIM(sn_rcv_wnum%clcat  ), ')'    
    298          WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%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  ), ')' 
    299304         WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 
    300305         WRITE(numout,*)'  sent fields (multiple ice categories)' 
     
    309314         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')' 
    310315         WRITE(numout,*)'      water level                     = ', TRIM(sn_snd_wlev%cldes  ), ' (', TRIM(sn_snd_wlev%clcat  ), ')'    
    311          WRITE(numout,*)'      mean sea level pressure         = ', TRIM(sn_rcv_mslp%cldes  ), ' (', TRIM(sn_rcv_mslp%clcat  ), ')'    
    312316         WRITE(numout,*)'      surface current to waves        = ', TRIM(sn_snd_crtw%cldes  ), ' (', TRIM(sn_snd_crtw%clcat  ), ')'    
    313317         WRITE(numout,*)'                      - referential   = ', sn_snd_crtw%clvref    
     
    525529      !                                                      ! Mean Sea Level Pressure   !    
    526530      !                                                      ! ------------------------- !    
    527       srcv(jpr_mslp)%clname = 'O_MSLP'     ;   IF( TRIM(sn_rcv_mslp%cldes  ) == 'coupled' )    srcv(jpr_mslp)%laction = .TRUE.    
     531      srcv(jpr_mslp)%clname = 'O_MSLP'  
     532      IF( TRIM(sn_rcv_mslp%cldes  ) == 'coupled' ) THEN  
     533         srcv(jpr_mslp)%laction = .TRUE.  
     534         cpl_mslp = .TRUE.  
     535      ENDIF 
    528536 
    529537      !                                                      ! ------------------------- ! 
     
    554562      ENDIF   
    555563      srcv(jpr_sdrftx)%clname = 'O_Sdrfx'    ! Stokes drift in the u direction   
    556       IF( TRIM(sn_rcv_sdrfx%cldes ) == 'coupled' )  THEN   
     564      srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction   
     565      IF( TRIM(sn_rcv_sdrft%cldes ) == 'coupled' )  THEN 
    557566         srcv(jpr_sdrftx)%laction = .TRUE.   
    558          cpl_sdrftx = .TRUE.   
    559       ENDIF   
    560       srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction   
    561       IF( TRIM(sn_rcv_sdrfy%cldes ) == 'coupled' )  THEN   
    562          srcv(jpr_sdrfty)%laction = .TRUE.   
    563          cpl_sdrfty = .TRUE.   
     567         cpl_sdrft = .TRUE. 
    564568      ENDIF   
    565569      srcv(jpr_wper)%clname = 'O_WPer'       ! mean wave period   
     
    568572         cpl_wper = .TRUE.   
    569573      ENDIF   
     574      srcv(jpr_wfreq)%clname = 'O_WFreq'     ! wave peak frequency   
     575      IF( TRIM(sn_rcv_wfreq%cldes ) == 'coupled' )  THEN   
     576         srcv(jpr_wfreq)%laction = .TRUE.   
     577         cpl_wfreq = .TRUE.   
     578      ENDIF 
    570579      srcv(jpr_wnum)%clname = 'O_WNum'       ! mean wave number   
    571580      IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' )  THEN   
     
    573582         cpl_wnum = .TRUE.   
    574583      ENDIF   
    575       srcv(jpr_wstrf)%clname = 'O_WStrf'     ! stress fraction adsorbed by the wave   
    576       IF( TRIM(sn_rcv_wstrf%cldes ) == 'coupled' )  THEN   
    577          srcv(jpr_wstrf)%laction = .TRUE.   
    578          cpl_wstrf = .TRUE.   
     584      srcv(jpr_tauoc)%clname = 'O_TauOce'     ! stress fraction adsorbed by the wave   
     585      IF( TRIM(sn_rcv_tauoc%cldes ) == 'coupled' )  THEN   
     586         srcv(jpr_tauoc)%laction = .TRUE.   
     587         cpl_tauoc = .TRUE.   
     588      ENDIF   
     589      srcv(jpr_tauwx)%clname = 'O_Tauwx'      ! ocean stress from wave in the x direction  
     590      srcv(jpr_tauwy)%clname = 'O_Tauwy'      ! ocean stress from wave in the y direction  
     591      IF( TRIM(sn_rcv_tauw%cldes ) == 'coupled' )  THEN   
     592         srcv(jpr_tauwx)%laction = .TRUE.   
     593         srcv(jpr_tauwy)%laction = .TRUE.   
     594         cpl_tauw = .TRUE. 
    579595      ENDIF   
    580596      srcv(jpr_wdrag)%clname = 'O_WDrag'     ! neutral surface drag coefficient   
     
    583599         cpl_wdrag = .TRUE.   
    584600      ENDIF   
     601      !  
     602      IF( srcv(jpr_tauoc)%laction .AND. srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction ) &  
     603            CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', &  
     604                                     '(sn_rcv_tauoc=coupled and sn_rcv_tauw=coupled)' )  
     605      ! 
    585606      !    
    586607      !                                                      ! ------------------------------- ! 
     
    964985      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case) 
    965986      !!---------------------------------------------------------------------- 
    966       USE zdf_oce,  ONLY : ln_zdfqiao 
     987      USE sbcflx ,  ONLY : ln_shelf_flx  
     988      USE sbcssm ,  ONLY : sbc_ssm_cpl  
     989      USE lib_fortran     ! distributed memory computing library 
    967990 
    968991      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index 
     
    10711094      ELSE                                                   !   No dynamical coupling   ! 
    10721095         !                                                   ! ========================= ! 
     1096         ! it is possible that the momentum is calculated from the winds (ln_shelf_flx) and a coupled drag coefficient  
     1097         IF( srcv(jpr_wdrag)%laction .AND. ln_shelf_flx .AND. ln_cdgw .AND. nn_drag == jp_std ) THEN  
     1098            DO jj = 1, jpjm1  
     1099               DO ji = 1, jpim1  
     1100                  ! here utau and vtau should contain the wind components as read from the forcing files,  
     1101                  ! and the wind module should already be properly calculated  
     1102                  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) ) * &  
     1103                                                                         utau(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) )  
     1104                  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) ) * &  
     1105                                                                         vtau(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) )  
     1106                  utau(ji,jj) = frcv(jpr_otx1)%z3(ji,jj,1)  
     1107                  vtau(ji,jj) = frcv(jpr_oty1)%z3(ji,jj,1)  
     1108               END DO  
     1109            END DO  
     1110            CALL lbc_lnk_multi( frcv(jpr_otx1)%z3(:,:,1), 'U', -1. , frcv(jpr_oty1)%z3(:,:,1), 'V', -1. , &  
     1111                                                             utau(:,:), 'U', -1. , vtau(:,:), 'V',  -1. )  
     1112            llnewtx = .TRUE.  
     1113         ELSE 
    10731114         frcv(jpr_otx1)%z3(:,:,1) = 0.e0                               ! here simply set to zero  
    10741115         frcv(jpr_oty1)%z3(:,:,1) = 0.e0                               ! an external read in a file can be added instead 
    10751116         llnewtx = .TRUE. 
     1117         ENDIF 
    10761118         ! 
    10771119      ENDIF 
     
    10931135            END DO 
    10941136            CALL lbc_lnk( frcv(jpr_taum)%z3(:,:,1), 'T', 1. ) 
     1137            IF( .NOT. srcv(jpr_otx1)%laction .AND. srcv(jpr_wdrag)%laction .AND. &  
     1138                                ln_shelf_flx .AND. ln_cdgw .AND. nn_drag == jp_std ) &  
     1139               taum(:,:) = frcv(jpr_taum)%z3(:,:,1 
    10951140            llnewtau = .TRUE. 
    10961141         ELSE 
     
    11071152      !                                                      ! ========================= ! 
    11081153      !                                                      !      10 m wind speed      !   (wndm) 
     1154      !                                                      !   include wave drag coef  !   (wndm) 
    11091155      !                                                      ! ========================= ! 
    11101156      ! 
     
    11171163!CDIR NOVERRCHK 
    11181164               DO ji = 1, jpi  
     1165                  IF( ln_shelf_flx ) THEN   ! the 10 wind module is properly calculated before if ln_shelf_flx  
     1166                     frcv(jpr_w10m)%z3(ji,jj,1) = wndm(ji,jj)  
     1167                  ELSE 
    11191168                  frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef ) 
     1169                  ENDIF 
    11201170               END DO 
    11211171            END DO 
     
    11271177      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN 
    11281178         ! 
     1179         ! if ln_wavcpl, the fields already contain the right information from forcing even if not ln_mixcpl 
    11291180         IF( ln_mixcpl ) THEN 
    1130             utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:) 
    1131             vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:) 
    1132             taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:) 
    1133             wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:) 
    1134          ELSE 
     1181            IF( srcv(jpr_otx1)%laction ) THEN  
     1182               utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:)  
     1183               vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:)  
     1184            ENDIF  
     1185            IF( srcv(jpr_taum)%laction )   &  
     1186               taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:)  
     1187            IF( srcv(jpr_w10m)%laction )   &  
     1188               wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:)  
     1189         ELSE IF( ll_purecpl ) THEN 
    11351190            utau(:,:) = frcv(jpr_otx1)%z3(:,:,1) 
    11361191            vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1) 
     
    11561211          IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields    
    11571212        
    1158           r1_grau = 1.e0 / (grav * rau0)               !* constant for optimization    
    1159           ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rpref ) * r1_grau    ! equivalent ssh (inverse barometer)    
     1213          !                                                  !* update the reference atmospheric pressure (if necessary)  
     1214          IF( ln_ref_apr )  rn_pref = glob_sum( frcv(jpr_mslp)%z3(:,:,1) * e1e2t(:,:) ) / tarea  
     1215         
     1216          ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rn_pref ) * r1_grau    ! equivalent ssh (inverse barometer)    
    11601217          apr   (:,:) =     frcv(jpr_mslp)%z3(:,:,1) !atmospheric pressure    
     1218          !  
     1219          CALL iom_put( "ssh_ib", ssh_ib )                                    !* output the inverse barometer ssh 
    11611220        
    1162           IF( kt == nit000 ) ssh_ibb(:,:) = ssh_ib(:,:)  ! correct this later (read from restart if possible)    
     1221          !                                         ! ---------------------------------------- !  
     1222          IF( kt == nit000 ) THEN                   !   set the forcing field at nit000 - 1    !  
     1223             !                                      ! ---------------------------------------- !  
     1224             !* Restart: read in restart file  
     1225             IF( ln_rstart .AND. iom_varid( numror, 'ssh_ibb', ldstop = .FALSE. ) > 0 ) THEN  
     1226                IF(lwp) WRITE(numout,*) 'sbc_cpl:   ssh_ibb read in the restart file'  
     1227                CALL iom_get( numror, jpdom_autoglo, 'ssh_ibb', ssh_ibb )   ! before inv. barometer ssh  
     1228             ELSE                                         !* no restart: set from nit000 values  
     1229                IF(lwp) WRITE(numout,*) 'sbc_cpl:   ssh_ibb set to nit000 values'  
     1230                ssh_ibb(:,:) = ssh_ib(:,:)  
     1231             ENDIF  
     1232          ENDIF  
     1233          !                                         ! ---------------------------------------- !  
     1234          IF( lrst_oce ) THEN                       !      Write in the ocean restart file     !  
     1235             !                                      ! ---------------------------------------- !  
     1236             IF(lwp) WRITE(numout,*)  
     1237             IF(lwp) WRITE(numout,*) 'sbc_cpl : ssh_ib written in ocean restart file at it= ', kt,' date= ', ndastp  
     1238             IF(lwp) WRITE(numout,*) '~~~~'  
     1239             CALL iom_rstput( kt, nitrst, numrow, 'ssh_ibb' , ssh_ib )  
     1240          ENDIF  
     1241         
     1242          ! Update mean ssh  
     1243          CALL sbc_ssm_cpl( kt ) 
    11631244      END IF    
    11641245      !   
    11651246      IF( ln_sdw ) THEN  ! Stokes Drift correction activated   
    11661247      !                                                      ! ========================= !    
    1167       !                                                      !       Stokes drift u      !   
     1248      !                                                      !       Stokes drift u,v    !   
    11681249      !                                                      ! ========================= !    
    1169          IF( srcv(jpr_sdrftx)%laction ) ut0sd(:,:) = frcv(jpr_sdrftx)%z3(:,:,1)   
    1170       !   
    1171       !                                                      ! ========================= !    
    1172       !                                                      !       Stokes drift v      !   
    1173       !                                                      ! ========================= !    
    1174          IF( srcv(jpr_sdrfty)%laction ) vt0sd(:,:) = frcv(jpr_sdrfty)%z3(:,:,1)   
     1250      IF( srcv(jpr_sdrftx)%laction .AND. srcv(jpr_sdrfty)%laction ) THEN  
     1251                                     ut0sd(:,:) = frcv(jpr_sdrftx)%z3(:,:,1)   
     1252                                     vt0sd(:,:) = frcv(jpr_sdrfty)%z3(:,:,1)   
     1253      ENDIF 
    11751254      !   
    11761255      !                                                      ! ========================= !    
     
    11851264      !   
    11861265      !                                                      ! ========================= !    
     1266      !                                                      !    Wave peak frequency    !   
     1267      !                                                      ! ========================= !    
     1268         IF( srcv(jpr_wfreq)%laction ) wfreq(:,:) = frcv(jpr_wfreq)%z3(:,:,1)   
     1269      !  
     1270      !                                                      ! ========================= ! 
    11871271      !                                                      !    Vertical mixing Qiao   !   
    11881272      !                                                      ! ========================= !    
     
    11901274        
    11911275         ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode   
    1192          IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction &   
    1193                                                                     .OR. srcv(jpr_hsig)%laction ) &   
     1276         IF( (srcv(jpr_sdrftx)%laction .AND. srcv(jpr_sdrfty)%laction) .OR. srcv(jpr_wper)%laction &   
     1277                                        .OR. srcv(jpr_hsig)%laction   .OR. srcv(jpr_wfreq)%laction) & 
    11941278            CALL sbc_stokes()   
    11951279      ENDIF   
     
    11971281      !                                                      ! Stress adsorbed by waves  !   
    11981282      !                                                      ! ========================= !    
    1199       IF( srcv(jpr_wstrf)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_wstrf)%z3(:,:,1)   
     1283      IF( srcv(jpr_tauoc)%laction .AND. ln_tauoc ) THEN  
     1284         tauoc_wave(:,:) = frcv(jpr_tauoc)%z3(:,:,1)  
     1285         ! cap the value of tauoc  
     1286         WHERE(tauoc_wave <   0.0 ) tauoc_wave = 1.0  
     1287         WHERE(tauoc_wave > 100.0 ) tauoc_wave = 1.0  
     1288      ENDIF  
     1289      !                                                      ! ========================= !    
     1290      !                                                      ! Stress component by waves !   
     1291      !                                                      ! ========================= !    
     1292      IF( srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction .AND. ln_tauw ) THEN  
     1293         tauw_x(:,:) = frcv(jpr_tauwx)%z3(:,:,1)  
     1294         tauw_y(:,:) = frcv(jpr_tauwy)%z3(:,:,1)  
     1295         ! cap the value of tauoc  
     1296         WHERE(tauw_x < -100.0 ) tauw_x = 0.0  
     1297         WHERE(tauw_x >  100.0 ) tauw_x = 0.0  
     1298         WHERE(tauw_y < -100.0 ) tauw_y = 0.0  
     1299         WHERE(tauw_y >  100.0 ) tauw_y = 0.0  
     1300      ENDIF 
    12001301        
    12011302      !                                                      ! ========================= !    
    1202       !                                                      ! Wave drag coefficient   !   
     1303      !                                                      !   Wave to ocean energy    ! 
    12031304      !                                                      ! ========================= !    
    1204       IF( srcv(jpr_wdrag)%laction .AND. ln_cdgw ) cdn_wave(:,:) = frcv(jpr_wdrag)%z3(:,:,1) 
     1305      IF( srcv(jpr_phioc)%laction .AND. ln_phioc ) THEN  
     1306         rn_crban(:,:) = 29.0 * frcv(jpr_phioc)%z3(:,:,1)  
     1307         WHERE( rn_crban <    0.0 ) rn_crban = 0.0  
     1308         WHERE( rn_crban > 1000.0 ) rn_crban = 1000.0  
     1309      ENDIF  
    12051310 
    12061311      !  Fields received by SAS when OASIS coupling 
     
    12741379               CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' ) 
    12751380            END SELECT 
    1276          ELSE 
     1381         ELSE IF( ll_purecpl ) THEN 
    12771382            zemp(:,:) = 0._wp 
    12781383         ENDIF 
     
    12821387         IF( srcv(jpr_cal)%laction )     zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1) 
    12831388          
    1284          IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:) 
    1285          ELSE                   ;   emp(:,:) =                              zemp(:,:) 
     1389         IF( ln_mixcpl .AND. ( srcv(jpr_oemp)%laction .OR. srcv(jpr_rain)%laction )) THEN  
     1390                                         emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:)  
     1391         ELSE IF( ll_purecpl ) THEN  ;   emp(:,:) = zemp(:,:) 
    12861392         ENDIF 
    12871393         ! 
     
    12981404            ENDIF 
    12991405         ENDIF 
    1300          IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:) 
    1301          ELSE                   ;   qns(:,:) =                              zqns(:,:) 
     1406         IF( ln_mixcpl .AND. ( srcv(jpr_qnsoce)%laction .OR. srcv(jpr_qnsmix)%laction )) THEN  
     1407                                          qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)  
     1408         ELSE IF( ll_purecpl ) THEN   ;   qns(:,:) = zqns(:,:) 
    13021409         ENDIF 
    13031410 
     
    13081415         ENDIF 
    13091416         IF( ln_dm2dc .AND. ln_cpl )   zqsr(:,:) = sbc_dcy( zqsr )   ! modify qsr to include the diurnal cycle 
    1310          IF( ln_mixcpl ) THEN   ;   qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:) 
    1311          ELSE                   ;   qsr(:,:) =                              zqsr(:,:) 
     1417         IF( ln_mixcpl .AND. ( srcv(jpr_qsroce)%laction .OR. srcv(jpr_qsrmix)%laction )) THEN  
     1418                                          qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:)  
     1419         ELSE IF( ll_purecpl ) THEN   ;   qsr(:,:) = zqsr(:,:) 
    13121420         ENDIF 
    13131421         ! 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90

    r8918 r10478  
    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 ) THEN 
    232288               wndm(ji,jj) = SQRT( zmod * zcoef ) 
     289               ENDIF 
    233290            END DO 
    234291         END DO 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r10473 r10478  
    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          &             ln_tauoc  , ln_stcor  , 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 
    134134         WRITE(numout,*) '              wave physics                               ln_wave     = ', ln_wave   
    135          WRITE(numout,*) '                 Stokes drift corr. to vert. velocity    ln_sdw      = ', ln_sdw   
    136          WRITE(numout,*) '                 wave modified ocean stress              ln_tauoc    = ', ln_tauoc   
    137          WRITE(numout,*) '                 Stokes coriolis term                    ln_stcor    = ', ln_stcor   
    138          WRITE(numout,*) '                 neutral drag coefficient (CORE, MFS)    ln_cdgw     = ', ln_cdgw 
    139135         WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis    = ', lk_oasis 
    140136         WRITE(numout,*) '              components of your executable              nn_components = ', nn_components 
    141137         WRITE(numout,*) '              Multicategory heat flux formulation (LIM3) nn_limflx   = ', nn_limflx 
     138         WRITE(numout,*) '              momentum formulation                       nn_drag     = ', nn_drag 
    142139         WRITE(numout,*) '           Misc. options of sbc : ' 
    143140         WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn 
     
    169166      IF ( nn_components == jp_iam_opa .AND. ln_cpl )   & 
    170167         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' ) 
    171       IF ( nn_components == jp_iam_opa .AND. ln_mixcpl )   & 
    172          &      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' ) 
    173170      IF ( ln_cpl .AND. .NOT. lk_oasis )    & 
    174171         &      CALL ctl_stop( 'STOP', 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' ) 
    175       IF( ln_mixcpl .AND. .NOT. lk_oasis )    & 
    176          &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires the cpp key key_oasis3' ) 
    177       IF( ln_mixcpl .AND. .NOT. ln_cpl )    & 
    178          &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires ln_cpl = T' ) 
    179       IF( ln_mixcpl .AND. nn_components /= jp_iam_nemo )    & 
    180          &      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' ) 
    181178 
    182179      !                              ! allocate sbc arrays 
     
    219216      IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk_core ) .AND. nn_components /= jp_iam_opa )   & 
    220217         &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' ) 
    221        
    222       IF ( ln_wave ) THEN 
    223          !Activated wave module but neither drag nor stokes drift activated   
    224          IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) )   THEN    
    225              CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F')   
    226          !drag coefficient read from wave model definable only with mfs bulk formulae and core 
    227          ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN        
    228              CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core') 
    229          ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN    
    230              CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 
    231          ENDIF 
    232       ELSE 
    233          IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) &    
    234             &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &   
    235             &                  'with drag coefficient (ln_cdgw =T) '  ,                        &   
    236             &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &   
    237             &                  'or ocean stress modification due to waves (ln_tauoc=T) ',      &   
    238             &                  'or Stokes-Coriolis term (ln_stcori=T)'  ) 
    239       ENDIF  
    240218      !                          ! Choice of the Surface Boudary Condition (set nsbc) 
    241       ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl 
     219      ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl .AND. .NOT. ln_wavcpl 
    242220      ! 
    243221      icpt = 0 
     
    271249         IF( nsbc == jp_mfs     )   WRITE(numout,*) '              MFS Bulk formulation' 
    272250         IF( nsbc == jp_none    )   WRITE(numout,*) '              OPA coupled to SAS via oasis' 
    273          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' 
    274253         IF( nn_components/= jp_iam_nemo )  & 
    275254            &                       WRITE(numout,*) '              + OASIS coupled SAS' 
     
    340319      !!              - updte the ice fraction : fr_i 
    341320      !!---------------------------------------------------------------------- 
     321      USE bdydta, ONLY: bdy_dta   
     322      ! 
    342323      INTEGER, INTENT(in) ::   kt       ! ocean time step 
    343324      !!--------------------------------------------------------------------- 
     
    360341      !                                            ! ---------------------------------------- ! 
    361342      ! 
    362       IF ( .NOT. lk_bdy ) then 
    363          IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc 
    364       ENDIF 
     343 
     344      IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc 
    365345                                                         ! (caution called before sbc_ssm) 
    366346      ! 
     
    397377      END SELECT 
    398378 
    399       IF( ln_mixcpl )        CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing 
    400  
    401       IF ( ln_wave .AND. ln_tauoc) THEN                 ! Wave stress subctracted   
    402             utau(:,:) = utau(:,:)*tauoc_wave(:,:)   
    403             vtau(:,:) = vtau(:,:)*tauoc_wave(:,:)   
    404             taum(:,:) = taum(:,:)*tauoc_wave(:,:)   
    405       !   
    406             SELECT CASE( nsbc )   
    407             CASE(  0,1,2,3,5,-1 )  ;   
    408                 IF(lwp .AND. kt == nit000 ) WRITE(numout,*) 'WARNING: You are subtracting the wave stress to the ocean. &   
    409                         & If not requested select ln_tauoc=.false'   
    410             END SELECT   
    411       !   
    412       END IF 
     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) 
    413384 
    414385      !                                            !==  Misc. Options  ==! 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90

    r8058 r10478  
    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_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r10473 r10478  
    1818   USE oce            ! ocean variables 
    1919   USE sbc_oce        ! Surface boundary condition: ocean fields 
    20    USE zdf_oce,  ONLY : ln_zdfqiao 
    2120   USE bdy_oce        ! open boundary condition variables 
    2221   USE domvvl         ! domain: variable volume layers 
     
    3332 
    3433   PUBLIC   sbc_stokes      ! routine called in sbccpl  
     34   PUBLIC   sbc_stress      ! routine called in sbcmod 
    3535   PUBLIC   sbc_wave        ! routine called in sbcmod  
    3636   PUBLIC   sbc_wave_init   ! routine called in sbcmod 
     
    3939   LOGICAL, PUBLIC ::   cpl_hsig   = .FALSE. 
    4040   LOGICAL, PUBLIC ::   cpl_phioc  = .FALSE. 
    41    LOGICAL, PUBLIC ::   cpl_sdrftx = .FALSE. 
    42    LOGICAL, PUBLIC ::   cpl_sdrfty = .FALSE. 
     41   LOGICAL, PUBLIC ::   cpl_sdrft  = .FALSE. 
    4342   LOGICAL, PUBLIC ::   cpl_wper   = .FALSE. 
     43   LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE. 
    4444   LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE. 
    45    LOGICAL, PUBLIC ::   cpl_wstrf  = .FALSE. 
     45   LOGICAL, PUBLIC ::   cpl_tauoc  = .FALSE. 
     46   LOGICAL, PUBLIC ::   cpl_tauw   = .FALSE. 
    4647   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 
    4753 
    4854   INTEGER ::   jpfld    ! number of files to read for stokes drift 
     
    5157   INTEGER ::   jp_hsw   ! index of significant wave hight      (m)      at T-point 
    5258   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 
    5360 
    5461   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
     
    5663   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao 
    5764   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  
    5867   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !: 
    5968   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 
    6071   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !: 
     72   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauw_x              !: 
     73   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauw_y              !: 
    6174   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !: 
    6275   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   div_sd              !: barotropic stokes drift divergence 
     
    6477   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd  , vsd  , wsd   !: Stokes drift velocities at u-, v- & w-points, resp. 
    6578 
    66    !! * Substitutions 
    6779#  include "vectopt_loop_substitute.h90" 
    6880   !!---------------------------------------------------------------------- 
     
    98110      CALL wrk_alloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd )  
    99111      ! 
    100       zfac = 2.0_wp * rpi / 16.0_wp 
    101       DO jj = 1, jpj               ! exp. wave number at t-point    (Eq. (19) in Breivick et al. (2014) ) 
    102          DO ji = 1, jpi 
    103             ! Stokes drift velocity estimated from Hs and Tmean 
    104             ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp ) 
    105             ! Stokes surface speed 
    106             tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj)) 
    107             ! Wavenumber scale 
    108             zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp ) 
    109          END DO 
    110       END DO 
    111       DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points 
    112          DO ji = 1, jpim1 
    113             zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) 
    114             zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) 
    115             ! 
    116             zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 
    117             zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 
    118          END DO 
    119       END DO 
     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 
    120146      ! 
    121147      !                       !==  horizontal Stokes Drift 3D velocity  ==! 
    122       DO jk = 1, jpkm1 
    123          DO jj = 2, jpjm1 
    124             DO ji = 2, jpim1 
    125                zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 
    126                zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 
    127                !                           
    128                zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth 
    129                zkh_v = zk_v(ji,jj) * zdep_v 
    130                !                                ! Depth attenuation 
    131                zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u ) 
    132                zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v ) 
    133                ! 
    134                usd(ji,jj,jk) = zda_u * zk_u(ji,jj) * umask(ji,jj,jk) 
    135                vsd(ji,jj,jk) = zda_v * zk_v(ji,jj) * vmask(ji,jj,jk) 
     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 
    136164            END DO 
    137165         END DO 
    138       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 
    139186      CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. ) 
    140187      ! 
     
    190237 
    191238 
     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 
    192276   SUBROUTINE sbc_wave( kt ) 
    193277      !!--------------------------------------------------------------------- 
     
    210294         CALL fld_read( kt, nn_fsbc, sf_cd )             ! read from external forcing 
    211295         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
    212       ENDIF 
    213  
    214       IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN    !==  Wave induced stress  ==! 
     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  ==! 
    215303         CALL fld_read( kt, nn_fsbc, sf_tauoc )          ! read wave norm stress from external forcing 
    216304         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 
    217       ENDIF 
    218  
    219       IF( ln_sdw )  THEN                           !==  Computation of the 3d Stokes Drift  ==!  
     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  ==!  
    220329         ! 
    221330         IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled 
    222331            CALL fld_read( kt, nn_fsbc, sf_sd )          ! read wave parameters from external forcing 
    223             IF( jp_hsw > 0 )   hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height 
    224             IF( jp_wmp > 0 )   wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
    225             IF( jp_usd > 0 )   ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
    226             IF( jp_vsd > 0 )   vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
    227          ENDIF 
    228          ! 
     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 
    229361         ! Read also wave number if needed, so that it is available in coupling routines 
    230362         IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN 
     
    235367         !                                         !==  Computation of the 3d Stokes Drift  ==!  
    236368         ! 
    237          IF( jpfld == 4 )   CALL sbc_stokes()            ! Calculate only if required fields are read 
    238          !                                               ! In coupled wave model-NEMO case the call is done after coupling 
     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 
    239374         ! 
    240375      ENDIF 
     
    261396      !! 
    262397      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files 
    263       TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i     ! array of namelist informations on the fields to read 
    264       TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  & 
    265                              &   sn_hsw, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read 
    266       ! 
    267       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc 
     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 
    268407      !!--------------------------------------------------------------------- 
    269408      ! 
     
    277416      IF(lwm) WRITE ( numond, namsbc_wave ) 
    278417      ! 
     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      ! 
    279478      IF( ln_cdgw ) THEN 
    280479         IF( .NOT. cpl_wdrag ) THEN 
    281480            ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    282             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
     481            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_cd structure' ) 
    283482            ! 
    284483                                   ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
    285484            IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
    286             CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' ) 
     485            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 
    287486         ENDIF 
    288487         ALLOCATE( cdn_wave(jpi,jpj) ) 
     
    290489 
    291490      IF( ln_tauoc ) THEN 
    292          IF( .NOT. cpl_wstrf ) THEN 
     491         IF( .NOT. cpl_tauoc ) THEN 
    293492            ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc 
    294             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
     493            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauoc structure' ) 
    295494            ! 
    296495                                    ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   ) 
    297496            IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 
    298             CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' ) 
     497            CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 
    299498         ENDIF 
    300499         ALLOCATE( tauoc_wave(jpi,jpj) ) 
    301500      ENDIF 
    302501 
    303       IF( ln_sdw ) THEN   ! Find out how many fields have to be read from file if not coupled 
    304          jpfld=0 
    305          jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0 
    306          IF( .NOT. cpl_sdrftx ) THEN 
     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 
    307537            jpfld  = jpfld + 1 
    308538            jp_usd = jpfld 
    309          ENDIF 
    310          IF( .NOT. cpl_sdrfty ) THEN 
    311539            jpfld  = jpfld + 1 
    312540            jp_vsd = jpfld 
    313541         ENDIF 
    314          IF( .NOT. cpl_hsig ) THEN 
     542         IF( .NOT. cpl_hsig .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 
    315543            jpfld  = jpfld + 1 
    316544            jp_hsw = jpfld 
    317545         ENDIF 
    318          IF( .NOT. cpl_wper ) THEN 
     546         IF( .NOT. cpl_wper .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 
    319547            jpfld  = jpfld + 1 
    320548            jp_wmp = jpfld 
    321549         ENDIF 
    322  
    323          ! Read from file only the non-coupled fields  
    324          IF( jpfld > 0 ) THEN 
    325             ALLOCATE( slf_i(jpfld) ) 
    326             IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd 
    327             IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd 
    328             IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw 
    329             IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp 
    330             ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift 
    331             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
    332             ! 
    333             DO ifpr= 1, jpfld 
    334                ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
    335                IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
    336             END DO 
    337             ! 
    338             CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' ) 
    339          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' ) 
     571         ! 
     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 
     576         ! 
     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 
    340581         ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) ) 
    341          ALLOCATE( hsw  (jpi,jpj)    , wmp  (jpi,jpj)     ) 
     582         ALLOCATE( wmp  (jpi,jpj)  ) 
     583         ALLOCATE( wfreq (jpi,jpj) ) 
    342584         ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     ) 
    343585         ALLOCATE( div_sd(jpi,jpj) ) 
     
    349591         IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN 
    350592            ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum 
    351             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wave structure' ) 
     593            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wn structure' ) 
    352594                                   ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
    353595            IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
    354             CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
     596            CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'read wave input', 'namsbc_wave' ) 
    355597         ENDIF 
    356598         ALLOCATE( wnum(jpi,jpj) ) 
     599      ENDIF 
     600 
     601      IF( ln_sdw .OR. ln_rough ) THEN 
     602         ALLOCATE( hsw  (jpi,jpj) ) 
    357603      ENDIF 
    358604      ! 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90

    r10473 r10478  
    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 
     
    3536   INTEGER , PUBLIC ::   nn_npc      !: non penetrative convective scheme call  frequency 
    3637   INTEGER , PUBLIC ::   nn_npcp     !: non penetrative convective scheme print frequency 
    37    LOGICAL , PUBLIC ::   ln_zdfqiao  !: Enhanced wave vertical mixing Qiao(2010) formulation flag 
    3838 
    3939 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r10473 r10478  
    2424   USE phycst         ! physical constants 
    2525   USE zdfmxl         ! mixed layer 
    26    USE sbcwave ,  ONLY: hsw   ! significant wave height   
     26   USE sbcwave 
    2727   ! 
    2828   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     
    4848   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points 
    4949   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 
    5055 
    5156   !                              !! ** Namelist  namzdf_gls  ** 
     
    6166   REAL(wp) ::   rn_emin           ! minimum value of TKE (m2/s2) 
    6267   REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 
    63    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 
    6469   REAL(wp) ::   rn_hsro           ! Minimum surface roughness 
    65    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)  
    6671 
    6772   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters 
     
    9297   REAL(wp) ::   rm7           =  0.0_wp 
    9398   REAL(wp) ::   rm8           =  0.318_wp 
    94    REAL(wp) ::   rtrans        =  0.1_wp 
     99   REAL(wp) ::   rtrans_default =  0.1_wp 
    95100   REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters 
    96    REAL(wp) ::   rsbc_tke1, rsbc_tke2, rfact_tke                  !     -           -           -        - 
    97    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                  !     -           -           -        - 
    98103   REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        - 
    99104   REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        - 
     
    117122      !!                ***  FUNCTION zdf_gls_alloc  *** 
    118123      !!---------------------------------------------------------------------- 
    119       ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     & 
    120          &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc ) 
     124      ALLOCATE( mxln(jpi,jpj,jpk) , zwall(jpi,jpj,jpk) ,     &  
     125         &      ustars2(jpi,jpj)  , ustarb2(jpi,jpj)   ,     &  
     126         &      rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) ,     &  
     127         &      rsbc_psi1(jpi,jpj), rtrans(jpi,jpj), STAT= zdf_gls_alloc ) 
    121128         ! 
    122129      IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc ) 
     
    163170      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp 
    164171 
     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 
    165185      IF( kt /= nit000 ) THEN   ! restore before value to compute tke 
    166186         avt (:,:,:) = avt_k (:,:,:) 
     
    200220         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 
    201221      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)   
    202          zhsro(:,:) = hsw(:,:) 
     222         zhsro = MAX(rn_frac_hs * hsw, rn_hsro) 
    203223      END SELECT 
    204224 
     
    314334      ! --------------------------------- 
    315335      ! 
     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 
    316401      SELECT CASE ( nn_bc_surf ) 
    317402      ! 
    318403      CASE ( 0 )             ! Dirichlet case 
    319404      ! First level 
    320       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) 
    321406      en(:,:,1) = MAX(en(:,:,1), rn_emin)  
    322407      z_elem_a(:,:,1) = en(:,:,1) 
     
    325410      !  
    326411      ! One level below 
    327       en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 
     412      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default * ((zhsro(:,:)+fsdepw(:,:,2)) & 
    328413          &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
    329414      en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
     
    334419      ! 
    335420      CASE ( 1 )             ! Neumann boundary condition on d(e)/dz 
    336       ! 
    337421      ! Dirichlet conditions at k=1 
    338       en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
     422      en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 
    339423      en(:,:,1)       = MAX(en(:,:,1), rn_emin)       
    340424      z_elem_a(:,:,1) = en(:,:,1) 
     
    346430      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    347431      z_elem_a(:,:,2) = 0._wp 
    348       zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 
     432      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:)) )) 
    349433      zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 
    350434           &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
     
    354438      ! 
    355439      END SELECT 
     440      ENDIF ! ln_phioc 
    356441 
    357442      ! Bottom boundary condition on tke 
     
    539624      ! 
    540625      CASE ( 0 )             ! Dirichlet boundary conditions 
    541       ! 
    542       ! Surface value 
    543       zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
    544       psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    545       z_elem_a(:,:,1) = psi(:,:,1) 
    546       z_elem_c(:,:,1) = 0._wp 
    547       z_elem_b(:,:,1) = 1._wp 
    548       ! 
    549       ! One level below 
    550       zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 
    551       zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 
    552       psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    553       z_elem_a(:,:,2) = 0._wp 
    554       z_elem_c(:,:,2) = 0._wp 
    555       z_elem_b(:,:,2) = 1._wp 
    556       !  
    557       ! 
     626         IF( ln_phioc ) THEN   ! Wave induced mixing case  
     627                               ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2  
     628                               ! balance between the production and the  
     629                               ! dissipation terms including the wave effect  
     630            IF( nn_wmix==jp_janssen ) THEN  
     631               ! Surface value  
     632               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic  
     633               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     634               z_elem_a(:,:,1) = psi(:,:,1)        
     635               z_elem_c(:,:,1) = 0._wp  
     636               z_elem_b(:,:,1) = 1._wp  
     637               !  
     638               ! One level below  
     639               zex1 = (rmm*ra_sf+rnn)  
     640               zex2 = (rmm*ra_sf)  
     641               zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2  
     642               psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1)  
     643               z_elem_a(:,:,2) = 0._wp  
     644               z_elem_c(:,:,2) = 0._wp  
     645               z_elem_b(:,:,2) = 1._wp  
     646               !   
     647               !  
     648            ELSE IF( nn_wmix==jp_craigbanner ) THEN  
     649               ! Surface value  
     650               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic  
     651               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     652               z_elem_a(:,:,1) = psi(:,:,1)  
     653               z_elem_c(:,:,1) = 0._wp  
     654               z_elem_b(:,:,1) = 1._wp  
     655               !  
     656               ! One level below  
     657               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdepw(:,:,2)/zhsro(:,:) )))  
     658               zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)  
     659               psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     660               z_elem_a(:,:,2) = 0._wp  
     661               z_elem_c(:,:,2) = 0._wp  
     662               z_elem_b(:,:,2) = 1._wp  
     663               !   
     664               !  
     665            ENDIF  
     666         ELSE                  ! Wave induced mixing case with default values  
     667            ! Surface value  
     668            zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic  
     669            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     670            z_elem_a(:,:,1) = psi(:,:,1)  
     671            z_elem_c(:,:,1) = 0._wp  
     672            z_elem_b(:,:,1) = 1._wp  
     673            !  
     674            ! One level below  
     675            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdepw(:,:,2)/zhsro(:,:) )))  
     676            zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)  
     677            psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     678            z_elem_a(:,:,2) = 0._wp  
     679            z_elem_c(:,:,2) = 0._wp  
     680            z_elem_b(:,:,2) = 1._wp  
     681            !   
     682            !  
     683         ENDIF 
    558684      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz 
    559       ! 
    560       ! Surface value: Dirichlet 
    561       zdep(:,:)       = zhsro(:,:) * rl_sf 
    562       psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    563       z_elem_a(:,:,1) = psi(:,:,1) 
    564       z_elem_c(:,:,1) = 0._wp 
    565       z_elem_b(:,:,1) = 1._wp 
    566       ! 
    567       ! Neumann condition at k=2 
    568       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    569       z_elem_a(:,:,2) = 0._wp 
    570       ! 
    571       ! Set psi vertical flux at the surface: 
    572       zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
    573       zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
    574       zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
    575       zdep(:,:) =  rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 
    576              & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 
    577       zflxs(:,:) = zdep(:,:) * zflxs(:,:) 
    578       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
    579  
    580       !    
    581       ! 
     685         IF( ln_phioc ) THEN  ! Wave induced mixing case with forced/coupled fields  
     686            IF( nn_wmix==jp_janssen ) THEN  
     687               !  
     688               zdep(:,:) = rl_sf * zhsro(:,:)  
     689               psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     690               z_elem_a(:,:,1) = psi(:,:,1)  
     691               z_elem_c(:,:,1) = 0._wp  
     692               z_elem_b(:,:,1) = 1._wp  
     693               !  
     694               ! Neumann condition at k=2  
     695               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b  
     696               z_elem_a(:,:,2) = 0._wp  
     697               !  
     698               ! Set psi vertical flux at the surface:  
     699               zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf)  
     700               zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) &  
     701                  &                   * en(:,:,1)**rmm * zdep  
     702               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)  
     703               !  
     704            ELSE IF( nn_wmix==jp_craigbanner ) THEN  
     705               ! Surface value: Dirichlet  
     706               zdep(:,:)       = zhsro(:,:) * rl_sf  
     707               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     708               z_elem_a(:,:,1) = psi(:,:,1)  
     709               z_elem_c(:,:,1) = 0._wp  
     710               z_elem_b(:,:,1) = 1._wp  
     711               !  
     712               ! Neumann condition at k=2  
     713               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b  
     714               z_elem_a(:,:,2) = 0._wp  
     715               !  
     716               ! Set psi vertical flux at the surface:  
     717               zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope  
     718               zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)  
     719               zflxs(:,:) = (rnn + rsbc_tke1(:,:) * (rnn + rmm*ra_sf) * zdep(:,:)) * &  
     720                          (1._wp + rsbc_tke1(:,:) * zdep(:,:))**(2._wp*rmm/3._wp-1_wp)  
     721               zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &  
     722                      & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.)  
     723               zflxs(:,:) = zdep(:,:) * zflxs(:,:)  
     724               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)  
     725               !     
     726            ENDIF  
     727         ELSE                 ! Wave induced mixing case with default values  
     728            ! Surface value: Dirichlet  
     729            zdep(:,:)       = zhsro(:,:) * rl_sf  
     730            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     731            z_elem_a(:,:,1) = psi(:,:,1)  
     732            z_elem_c(:,:,1) = 0._wp  
     733            z_elem_b(:,:,1) = 1._wp  
     734            !  
     735            ! Neumann condition at k=2  
     736            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b  
     737            z_elem_a(:,:,2) = 0._wp  
     738            !  
     739            ! Set psi vertical flux at the surface:  
     740            zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope  
     741            zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)  
     742            zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * &  
     743                       (1._wp + rsbc_tke1_default * zdep(:,:))**(2._wp*rmm/3._wp-1_wp)  
     744            zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &  
     745                   & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.)  
     746            zflxs(:,:) = zdep(:,:) * zflxs(:,:)  
     747            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)  
     748            !     
     749            !  
     750         ENDIF 
    582751      END SELECT 
    583752 
     
    8701039      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 
    8711040         &            rn_clim_galp, ln_sigpsi, rn_hsro,      & 
    872          &            rn_crban, rn_charn, rn_frac_hs,        & 
     1041         &            rn_crban_default, rn_charn, rn_frac_hs,& 
    8731042         &            nn_bc_surf, nn_bc_bot, nn_z0_met,      & 
    8741043         &            nn_stab_func, nn_clos 
     
    8981067         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot 
    8991068         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi 
    900          WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban 
     1069         WRITE(numout,*) '      Craig and Banner coefficient (default)        rn_crban       = ', rn_crban_default 
    9011070         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn 
    9021071         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met 
     
    9151084      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' )   
    9161085      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' )   
    917       IF( nn_z0_met == 3 .AND. .NOT.ln_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )   
     1086      IF( nn_z0_met == 3 .AND. .NOT.ln_rough ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_rough=T' )   
     1087      IF( nn_z0_met .NE. 3 .AND. ln_rough ) THEN  
     1088         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')   
     1089         nn_z0_met = 3   
     1090      ENDIF 
    9181091      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' )   
    9191092      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 
     
    10891262               &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 
    10901263 
    1091       IF ( rn_crban==0._wp ) THEN 
     1264      IF( .NOT. ln_phioc .AND. rn_crban_default==0._wp ) THEN 
    10921265         rl_sf = vkarmn 
    10931266      ELSE 
     
    11281301      rc03  = rc02 * rc0 
    11291302      rc04  = rc03 * rc0 
    1130       rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking 
    1131       rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking  
    1132       zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
    1133       rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer  
     1303      rsbc_tke1_default = -3._wp/2._wp*rn_crban_default*ra_sf*rl_sf  
     1304      rsbc_tke2         = rdt * rn_crban_default / rl_sf  
     1305      zcr               = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp )  
     1306      rtrans_default    = 0.2_wp / zcr 
    11341307      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness 
    11351308      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness  
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90

    r10473 r10478  
    5454      !! 
    5555      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,       &   
    57          &        ln_zdfqiao 
     56         &        ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp 
    5857      !!---------------------------------------------------------------------- 
    5958 
     
    8483         WRITE(numout,*) '      npc call  frequency                 nn_npc    = ', nn_npc 
    8584         WRITE(numout,*) '      npc print frequency                 nn_npcp   = ', nn_npcp 
    86          WRITE(numout,*) '      Qiao formulation flag               ln_zdfqiao=', ln_zdfqiao 
    8785      ENDIF 
    8886 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/step.F90

    r10473 r10478  
    106106      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    107107      IF( lk_tide    )   CALL sbc_tide( kstp ) 
    108       IF( lk_bdy     )  THEN 
    109          IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib  
    110                          CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries 
    111       ENDIF 
    112                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice) 
     108                         CALL sbc     ( kstp )        ! Sea Boundary Condition (including sea-ice) 
    113109                                                      ! clem: moved here for bdy ice purpose 
    114110      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
Note: See TracChangeset for help on using the changeset viewer.