Changeset 7471


Ignore:
Timestamp:
2016-12-07T13:15:34+01:00 (4 years ago)
Author:
jcastill
Message:

Version as merged to the trunk during the Nov-2016 merge party, equivalent to branches/UKMO/r5936_INGV1_WAVE-coupling@7360

Location:
branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM
Files:
2 added
14 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref

    r5578 r7471  
    264264                           !     =1 global mean of e-p-r set to zero at each time step 
    265265                           !     =2 annual global mean of e-p-r set to zero 
    266    ln_wave = .false.       !  Activate coupling with wave (either Stokes Drift or Drag coefficient, or both)  (T => fill namsbc_wave) 
    267    ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave) 
    268    ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave) 
     266   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)   
    269271   nn_lsm  = 0             !  =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) , 
    270272                           !  =1:n number of iterations of land/sea mask application for input fields (fill land/sea mask filename field) 
     
    363365   sn_snd_crt    =       'none'                 ,    'no'    , 'spherical' , 'eastward-northward' ,  'T' 
    364366   sn_snd_co2    =       'coupled'              ,    'no'    ,     ''      ,         ''           ,   '' 
     367   sn_snd_crtw   =       'none'                 ,    'no'    ,     ''      ,         ''           , 'U,V'  
     368   sn_snd_ifrac  =       'none'                 ,    'no'    ,     ''      ,         ''           ,   ''  
     369   sn_snd_wlev   =       'none'                 ,    'no'    ,     ''      ,         ''           ,   ''  
    365370! receive 
    366371   sn_rcv_w10m   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
     
    374379   sn_rcv_cal    =       'coupled'              ,    'no'    ,     ''      ,         ''          ,   '' 
    375380   sn_rcv_co2    =       'coupled'              ,    'no'    ,     ''      ,         ''          ,   '' 
     381   sn_rcv_hsig   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
     382   sn_rcv_iceflx =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
     383   sn_rcv_mslp   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
     384   sn_rcv_phioc  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
     385   sn_rcv_sdrfx  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
     386   sn_rcv_sdrfy  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
     387   sn_rcv_wper   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
     388   sn_rcv_wnum   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
     389   sn_rcv_wstrf  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
     390   sn_rcv_wdrag  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
    376391! 
    377392   nn_cplmodel   =     1     !  Maximum number of models to/from which NEMO is potentialy sending/receiving data 
     
    921936   ln_zdfexp   = .false.   !  time-stepping: split-explicit (T) or implicit (F) time stepping 
    922937   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) 
    923939/ 
    924940!----------------------------------------------------------------------- 
     
    12701286!              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
    12711287!              !             !  (if <0  months)  !   name       !   (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      ! 
    1272    sn_cdg      =  'cdg_wave' ,        1          , 'drag_coeff' ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
     1288   sn_cdg      =  'sdw_wave' ,        1          , 'drag_coeff' ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    12731289   sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    12741290   sn_vsd      =  'sdw_wave' ,        1          , 'v_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    1275    sn_wn       =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    1276 ! 
    1277    cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files 
     1291   sn_swh      =  'sdw_wave' ,        1          , 'hs'         ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
     1292   sn_wmp      =  'sdw_wave' ,        1          , 'wmp'        ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
     1293   sn_wnum     =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
     1294   sn_tauoc    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
     1295!  
     1296   cn_dir  = './'  !  root directory for the location of drag coefficient files 
    12781297/ 
    12791298!----------------------------------------------------------------------- 
  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/CONFIG/cfg.txt

    r7470 r7471  
    1111ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
    1212ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 
     13GYRE_LONG OPA_SRC 
  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90

    r7470 r7471  
    6565   INTEGER                    ::   nsnd         ! total number of fields sent  
    6666   INTEGER                    ::   ncplmodel    ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 
    67    INTEGER, PUBLIC, PARAMETER ::   nmaxfld=50   ! Maximum number of coupling fields 
     67   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=55   ! Maximum number of coupling fields 
    6868   INTEGER, PUBLIC, PARAMETER ::   nmaxcat=5    ! Maximum number of coupling fields 
    6969   INTEGER, PUBLIC, PARAMETER ::   nmaxcpl=5    ! Maximum number of coupling fields 
  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r7470 r7471  
    6565   LOGICAL , PUBLIC ::   ln_cdgw        !: true if neutral drag coefficient from wave model 
    6666   LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model 
     67   LOGICAL , PUBLIC ::   ln_tauoc       !: true if normalized stress from wave is used  
     68   LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used  
    6769   ! 
    6870   LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs 
  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r7470 r7471  
    737737 
    738738      !! Neutral coefficients at 10m: 
    739       IF( ln_cdgw ) THEN      ! wave drag case 
     739      IF( ln_wave .AND. ln_cdgw ) THEN      ! wave drag case 
    740740         cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) 
    741741         ztmp0   (:,:) = cdn_wave(:,:) 
     
    783783         END IF 
    784784        
    785          IF( ln_cdgw ) THEN      ! surface wave case 
     785         IF( ln_wave .AND. ln_cdgw ) THEN      ! surface wave case 
    786786            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u )  
    787787            Cd      = sqrt_Cd * sqrt_Cd 
  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r7470 r7471  
    2323   USE sbcapr 
    2424   USE sbcdcy          ! surface boundary condition: diurnal cycle 
     25   USE sbcwave         ! surface boundary condition: waves 
    2526   USE phycst          ! physical constants 
    2627#if defined key_lim3 
     
    105106   INTEGER, PARAMETER ::   jpr_e3t1st = 41            ! first T level thickness  
    106107   INTEGER, PARAMETER ::   jpr_fraqsr = 42            ! fraction of solar net radiation absorbed in the first ocean level 
    107    INTEGER, PARAMETER ::   jprcv      = 42            ! total number of fields received 
     108   INTEGER, PARAMETER ::   jpr_mslp   = 43            ! mean sea level pressure   
     109   INTEGER, PARAMETER ::   jpr_hsig   = 44            ! Hsig   
     110   INTEGER, PARAMETER ::   jpr_phioc  = 45            ! Wave=>ocean energy flux   
     111   INTEGER, PARAMETER ::   jpr_sdrftx = 46            ! Stokes drift on grid 1   
     112   INTEGER, PARAMETER ::   jpr_sdrfty = 47            ! Stokes drift on grid 2   
     113   INTEGER, PARAMETER ::   jpr_wper   = 48            ! Mean wave period  
     114   INTEGER, PARAMETER ::   jpr_wnum   = 49            ! Mean wavenumber  
     115   INTEGER, PARAMETER ::   jpr_wstrf  = 50            ! Stress fraction adsorbed by waves  
     116   INTEGER, PARAMETER ::   jpr_wdrag  = 51            ! Neutral surface drag coefficient  
     117   INTEGER, PARAMETER ::   jprcv      = 51            ! total number of fields received 
    108118 
    109119   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere 
     
    135145   INTEGER, PARAMETER ::   jps_e3t1st = 27            ! first level depth (vvl) 
    136146   INTEGER, PARAMETER ::   jps_fraqsr = 28            ! fraction of solar net radiation absorbed in the first ocean level 
    137    INTEGER, PARAMETER ::   jpsnd      = 28            ! total number of fields sended 
     147   INTEGER, PARAMETER ::   jps_ficet  = 29            ! total ice fraction    
     148   INTEGER, PARAMETER ::   jps_ocxw   = 30            ! currents on grid 1    
     149   INTEGER, PARAMETER ::   jps_ocyw   = 31            ! currents on grid 2  
     150   INTEGER, PARAMETER ::   jps_wlev   = 32            ! water level   
     151   INTEGER, PARAMETER ::   jpsnd      = 32            ! total number of fields sent 
    138152 
    139153   !                                                         !!** namelist namsbc_cpl ** 
     
    149163   ! Received from the atmosphere                     ! 
    150164   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf 
    151    TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                         
     165   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp                             
     166   ! Send to waves   
     167   TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev   
     168   ! Received from waves   
     169   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 
    152170   ! Other namelist parameters                        ! 
    153171   INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 
     
    161179 
    162180   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky) 
     181    
     182   REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2]   
     183   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0 
    163184 
    164185   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument 
     
    179200      !!             ***  FUNCTION sbc_cpl_alloc  *** 
    180201      !!---------------------------------------------------------------------- 
    181       INTEGER :: ierr(3) 
     202      INTEGER :: ierr(4) 
    182203      !!---------------------------------------------------------------------- 
    183204      ierr(:) = 0 
     
    190211      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) ) 
    191212      ! 
     213      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) )   
     214   
    192215      sbc_cpl_alloc = MAXVAL( ierr ) 
    193216      IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc ) 
     
    216239      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos 
    217240      !! 
    218       NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,      & 
    219          &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,      & 
    220          &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx,   & 
    221          &                  sn_rcv_co2 , nn_cplmodel  , ln_usecplmask 
     241      NAMELIST/namsbc_cpl/  sn_snd_temp , sn_snd_alb  , sn_snd_thick , sn_snd_crt   , sn_snd_co2,      &   
     242         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,      &   
     243         &                  sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev  , sn_rcv_hsig  , sn_rcv_phioc ,   &   
     244         &                  sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_wstrf ,   &  
     245         &                  sn_rcv_wdrag, sn_rcv_qns  , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   ,   &  
     246         &                  sn_rcv_iceflx,sn_rcv_co2  , nn_cplmodel  , ln_usecplmask, sn_rcv_mslp 
    222247      !!--------------------------------------------------------------------- 
    223248      ! 
     
    260285         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 
    261286         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')' 
     287         WRITE(numout,*)'      significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')'   
     288         WRITE(numout,*)'      wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')'   
     289         WRITE(numout,*)'      Surface Stokes drift grid u     = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')'   
     290         WRITE(numout,*)'      Surface Stokes drift grid v     = ', TRIM(sn_rcv_sdrfy%cldes ), ' (', TRIM(sn_rcv_sdrfy%clcat ), ')'   
     291         WRITE(numout,*)'      Mean wave period                = ', TRIM(sn_rcv_wper%cldes  ), ' (', TRIM(sn_rcv_wper%clcat  ), ')'   
     292         WRITE(numout,*)'      Mean wave number                = ', TRIM(sn_rcv_wnum%cldes  ), ' (', TRIM(sn_rcv_wnum%clcat  ), ')'   
     293         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')'   
     294         WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 
    262295         WRITE(numout,*)'  sent fields (multiple ice categories)' 
    263296         WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')' 
    264297         WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')' 
    265298         WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')' 
     299         WRITE(numout,*)'      total ice fraction              = ', TRIM(sn_snd_ifrac%cldes ), ' (', TRIM(sn_snd_ifrac%clcat ), ')' 
    266300         WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')' 
    267301         WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref  
     
    269303         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd 
    270304         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')' 
     305         WRITE(numout,*)'      water level                     = ', TRIM(sn_snd_wlev%cldes  ), ' (', TRIM(sn_snd_wlev%clcat  ), ')'   
     306         WRITE(numout,*)'      mean sea level pressure         = ', TRIM(sn_rcv_mslp%cldes  ), ' (', TRIM(sn_rcv_mslp%clcat  ), ')'   
     307         WRITE(numout,*)'      surface current to waves        = ', TRIM(sn_snd_crtw%cldes  ), ' (', TRIM(sn_snd_crtw%clcat  ), ')'   
     308         WRITE(numout,*)'                      - referential   = ', sn_snd_crtw%clvref   
     309         WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor   
     310         WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd 
    271311         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel 
    272312         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask 
     
    307347      !  
    308348      ! Vectors: change of sign at north fold ONLY if on the local grid 
     349      IF( TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM(sn_rcv_tau%cldes ) == 'oce and ice') THEN ! avoid working with the atmospheric fields if they are not coupled 
    309350      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1. 
    310351       
     
    374415         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp. 
    375416      ENDIF 
     417      ENDIF 
    376418        
    377419      !                                                      ! ------------------------- ! 
     
    470512      !                                                      ! ------------------------- ! 
    471513      srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE. 
     514       
     515      !                                                      ! ------------------------- !   
     516      !                                                      ! Mean Sea Level Pressure   !   
     517      !                                                      ! ------------------------- !   
     518      srcv(jpr_mslp)%clname = 'O_MSLP'     ;   IF( TRIM(sn_rcv_mslp%cldes  ) == 'coupled' )    srcv(jpr_mslp)%laction = .TRUE.   
     519       
    472520      !                                                      ! ------------------------- ! 
    473521      !                                                      !   topmelt and botmelt     !    
     
    483531         srcv(jpr_topm:jpr_botm)%laction = .TRUE. 
    484532      ENDIF 
     533      !                                                      ! ------------------------- !  
     534      !                                                      !      Wave breaking        !      
     535      !                                                      ! ------------------------- !   
     536      srcv(jpr_hsig)%clname  = 'O_Hsigwa'    ! significant wave height  
     537      IF( TRIM(sn_rcv_hsig%cldes  ) == 'coupled' )  THEN  
     538         srcv(jpr_hsig)%laction = .TRUE.  
     539         cpl_hsig = .TRUE.  
     540      ENDIF  
     541      srcv(jpr_phioc)%clname = 'O_PhiOce'    ! wave to ocean energy  
     542      IF( TRIM(sn_rcv_phioc%cldes ) == 'coupled' )  THEN  
     543         srcv(jpr_phioc)%laction = .TRUE.  
     544         cpl_phioc = .TRUE.  
     545      ENDIF  
     546      srcv(jpr_sdrftx)%clname = 'O_Sdrfx'    ! Stokes drift in the u direction  
     547      IF( TRIM(sn_rcv_sdrfx%cldes ) == 'coupled' )  THEN  
     548         srcv(jpr_sdrftx)%laction = .TRUE.  
     549         cpl_sdrftx = .TRUE.  
     550      ENDIF  
     551      srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction  
     552      IF( TRIM(sn_rcv_sdrfy%cldes ) == 'coupled' )  THEN  
     553         srcv(jpr_sdrfty)%laction = .TRUE.  
     554         cpl_sdrfty = .TRUE.  
     555      ENDIF  
     556      srcv(jpr_wper)%clname = 'O_WPer'       ! mean wave period  
     557      IF( TRIM(sn_rcv_wper%cldes  ) == 'coupled' )  THEN  
     558         srcv(jpr_wper)%laction = .TRUE.  
     559         cpl_wper = .TRUE.  
     560      ENDIF  
     561      srcv(jpr_wnum)%clname = 'O_WNum'       ! mean wave number  
     562      IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' )  THEN  
     563         srcv(jpr_wnum)%laction = .TRUE.  
     564         cpl_wnum = .TRUE.  
     565      ENDIF  
     566      srcv(jpr_wstrf)%clname = 'O_WStrf'     ! stress fraction adsorbed by the wave  
     567      IF( TRIM(sn_rcv_wstrf%cldes ) == 'coupled' )  THEN  
     568         srcv(jpr_wstrf)%laction = .TRUE.  
     569         cpl_wstrf = .TRUE.  
     570      ENDIF  
     571      srcv(jpr_wdrag)%clname = 'O_WDrag'     ! neutral surface drag coefficient  
     572      IF( TRIM(sn_rcv_wdrag%cldes ) == 'coupled' )  THEN  
     573         srcv(jpr_wdrag)%laction = .TRUE.  
     574         cpl_wdrag = .TRUE.  
     575      ENDIF  
     576      !   
    485577      !                                                      ! ------------------------------- ! 
    486578      !                                                      !   OPA-SAS coupling - rcv by opa !    
     
    637729      !                                                      ! ------------------------- ! 
    638730      ssnd(jps_fice)%clname = 'OIceFrc' 
     731      ssnd(jps_ficet)%clname = 'OIceFrcT' 
    639732      ssnd(jps_hice)%clname = 'OIceTck' 
    640733      ssnd(jps_hsnw)%clname = 'OSnwTck' 
     
    645738      ENDIF 
    646739       
     740      IF (TRIM( sn_snd_ifrac%cldes )  == 'coupled') ssnd(jps_ficet)%laction = .TRUE. 
     741 
    647742      SELECT CASE ( TRIM( sn_snd_thick%cldes ) ) 
    648743      CASE( 'none'         )       ! nothing to do 
     
    665760      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1' 
    666761      ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1' 
     762      ssnd(jps_ocxw)%clname = 'O_OCurxw'   
     763      ssnd(jps_ocyw)%clname = 'O_OCuryw' 
    667764      ! 
    668765      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold 
     
    685782      END SELECT 
    686783 
     784      ssnd(jps_ocxw:jps_ocyw)%nsgn = -1.   ! vectors: change of the sign at the north fold   
     785              
     786      IF( sn_snd_crtw%clvgrd == 'U,V' ) THEN   
     787         ssnd(jps_ocxw)%clgrid = 'U' ; ssnd(jps_ocyw)%clgrid = 'V'   
     788      ELSE IF( sn_snd_crtw%clvgrd /= 'T' ) THEN   
     789         CALL ctl_stop( 'sn_snd_crtw%clvgrd must be equal to T' )   
     790      ENDIF   
     791      IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) ssnd(jps_ocxw:jps_ocyw)%nsgn = 1.   
     792      SELECT CASE( TRIM( sn_snd_crtw%cldes ) )   
     793         CASE( 'none'                 )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .FALSE.   
     794         CASE( 'oce only'             )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .TRUE.   
     795         CASE( 'weighted oce and ice' )   !   nothing to do   
     796         CASE( 'mixed oce-ice'        )   ; ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.   
     797         CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crtw%cldes' )   
     798      END SELECT   
     799       
    687800      !                                                      ! ------------------------- ! 
    688801      !                                                      !          CO2 flux         ! 
    689802      !                                                      ! ------------------------- ! 
    690803      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE. 
     804 
     805      !                                                      ! ------------------------- !   
     806      !                                                      !     Sea surface height    !   
     807      !                                                      ! ------------------------- !   
     808      ssnd(jps_wlev)%clname = 'O_Wlevel' ;  IF( TRIM(sn_snd_wlev%cldes) == 'coupled' )   ssnd(jps_wlev)%laction = .TRUE.  
    691809 
    692810      !                                                      ! ------------------------------- ! 
     
    783901      IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   & 
    784902         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' ) 
    785       ncpl_qsr_freq = 86400 / ncpl_qsr_freq 
     903      IF( ln_dm2dc .AND. ln_cpl ) ncpl_qsr_freq = 86400 / ncpl_qsr_freq 
    786904 
    787905      CALL wrk_dealloc( jpi,jpj, zacs, zaos ) 
     
    837955      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case) 
    838956      !!---------------------------------------------------------------------- 
     957      USE zdf_oce,  ONLY : ln_zdfqiao 
     958 
    839959      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index 
    840960      INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation  
     
    9961116      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1) 
    9971117#endif 
     1118      !   
     1119      !                                                      ! ========================= !   
     1120      !                                                      ! Mean Sea Level Pressure   !   (taum)   
     1121      !                                                      ! ========================= !   
     1122      !   
     1123      IF( srcv(jpr_mslp)%laction ) THEN                    ! UKMO SHELF effect of atmospheric pressure on SSH   
     1124          IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields   
     1125       
     1126          r1_grau = 1.e0 / (grav * rau0)               !* constant for optimization   
     1127          ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rpref ) * r1_grau    ! equivalent ssh (inverse barometer)   
     1128          apr   (:,:) =     frcv(jpr_mslp)%z3(:,:,1) !atmospheric pressure   
     1129       
     1130          IF( kt == nit000 ) ssh_ibb(:,:) = ssh_ib(:,:)  ! correct this later (read from restart if possible)   
     1131      END IF   
     1132      !  
     1133      IF( ln_sdw ) THEN  ! Stokes Drift correction activated  
     1134      !                                                      ! ========================= !   
     1135      !                                                      !       Stokes drift u      !  
     1136      !                                                      ! ========================= !   
     1137         IF( srcv(jpr_sdrftx)%laction ) zusd2dt(:,:) = frcv(jpr_sdrftx)%z3(:,:,1)  
     1138      !  
     1139      !                                                      ! ========================= !   
     1140      !                                                      !       Stokes drift v      !  
     1141      !                                                      ! ========================= !   
     1142         IF( srcv(jpr_sdrfty)%laction ) zvsd2dt(:,:) = frcv(jpr_sdrfty)%z3(:,:,1)  
     1143      !  
     1144      !                                                      ! ========================= !   
     1145      !                                                      !      Wave mean period     !  
     1146      !                                                      ! ========================= !   
     1147         IF( srcv(jpr_wper)%laction ) wmp(:,:) = frcv(jpr_wper)%z3(:,:,1)  
     1148      !  
     1149      !                                                      ! ========================= !   
     1150      !                                                      !  Significant wave height  !  
     1151      !                                                      ! ========================= !   
     1152         IF( srcv(jpr_hsig)%laction ) swh(:,:) = frcv(jpr_hsig)%z3(:,:,1)  
     1153      !  
     1154      !                                                      ! ========================= !   
     1155      !                                                      !    Vertical mixing Qiao   !  
     1156      !                                                      ! ========================= !   
     1157         IF( srcv(jpr_wnum)%laction .AND. ln_zdfqiao ) wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1)  
     1158       
     1159         ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode  
     1160         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction &  
     1161                                                                    .OR. srcv(jpr_hsig)%laction ) THEN  
     1162            CALL sbc_stokes()  
     1163            IF( ln_zdfqiao .AND. .NOT. srcv(jpr_wnum)%laction ) CALL sbc_qiao()  
     1164         ENDIF  
     1165         IF( ln_zdfqiao .AND. srcv(jpr_wnum)%laction ) CALL sbc_qiao()  
     1166      ENDIF  
     1167      !                                                      ! ========================= !   
     1168      !                                                      ! Stress adsorbed by waves  !  
     1169      !                                                      ! ========================= !   
     1170      IF( srcv(jpr_wstrf)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_wstrf)%z3(:,:,1)  
     1171       
     1172      !                                                      ! ========================= !   
     1173      !                                                      !   Wave drag coefficient   !  
     1174      !                                                      ! ========================= !   
     1175      IF( srcv(jpr_wdrag)%laction .AND. ln_cdgw ) cdn_wave(:,:) = frcv(jpr_wdrag)%z3(:,:,1) 
    9981176 
    9991177      !  Fields received by SAS when OASIS coupling 
     
    20192197      ENDIF 
    20202198      ! 
     2199      !                                                      ! ------------------------- !   
     2200      !                                                      !  Surface current to waves !   
     2201      !                                                      ! ------------------------- !   
     2202      IF( ssnd(jps_ocxw)%laction .OR. ssnd(jps_ocyw)%laction ) THEN   
     2203          !       
     2204          !                                                  j+1  j     -----V---F   
     2205          ! surface velocity always sent from T point                    !       |   
     2206          !                                                       j      |   T   U   
     2207          !                                                              |       |   
     2208          !                                                   j   j-1   -I-------|   
     2209          !                                               (for I)        |       |   
     2210          !                                                             i-1  i   i   
     2211          !                                                              i      i+1 (for I)   
     2212          SELECT CASE( TRIM( sn_snd_crtw%cldes ) )   
     2213          CASE( 'oce only'             )      ! C-grid ==> T   
     2214             DO jj = 2, jpjm1   
     2215                DO ji = fs_2, fs_jpim1   ! vector opt.   
     2216                   zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )   
     2217                   zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji , jj-1,1) )    
     2218                END DO   
     2219             END DO   
     2220          CASE( 'weighted oce and ice' )      
     2221             SELECT CASE ( cp_ice_msh )   
     2222             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T   
     2223                DO jj = 2, jpjm1   
     2224                   DO ji = fs_2, fs_jpim1   ! vector opt.   
     2225                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un (ji-1,jj  ,1) ) * zfr_l(ji,jj)     
     2226                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn (ji  ,jj-1,1) ) * zfr_l(ji,jj)   
     2227                      zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)   
     2228                      zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)   
     2229                   END DO   
     2230                END DO   
     2231             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T   
     2232                DO jj = 2, jpjm1   
     2233                   DO ji = 2, jpim1   ! NO vector opt.   
     2234                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)     
     2235                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)     
     2236                      zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &   
     2237                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)   
     2238                      zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &   
     2239                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)   
     2240                   END DO   
     2241                END DO   
     2242             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T   
     2243                DO jj = 2, jpjm1   
     2244                   DO ji = 2, jpim1   ! NO vector opt.   
     2245                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)     
     2246                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)     
     2247                      zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &   
     2248                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)   
     2249                      zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &   
     2250                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)   
     2251                   END DO   
     2252                END DO   
     2253             END SELECT   
     2254             CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )   
     2255          CASE( 'mixed oce-ice'        )   
     2256             SELECT CASE ( cp_ice_msh )   
     2257             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T   
     2258                DO jj = 2, jpjm1   
     2259                   DO ji = fs_2, fs_jpim1   ! vector opt.   
     2260                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
     2261                         &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)   
     2262                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   &   
     2263                         &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)   
     2264                   END DO   
     2265                END DO   
     2266             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T   
     2267                DO jj = 2, jpjm1   
     2268                   DO ji = 2, jpim1   ! NO vector opt.   
     2269                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &      
     2270                         &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &   
     2271                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)   
     2272                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   &    
     2273                         &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &   
     2274                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)   
     2275                   END DO   
     2276                END DO   
     2277             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T   
     2278                DO jj = 2, jpjm1   
     2279                   DO ji = 2, jpim1   ! NO vector opt.   
     2280                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &      
     2281                         &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &   
     2282                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)   
     2283                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   &    
     2284                         &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &   
     2285                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)   
     2286                   END DO   
     2287                END DO   
     2288             END SELECT   
     2289          END SELECT   
     2290         CALL lbc_lnk( zotx1, ssnd(jps_ocxw)%clgrid, -1. )   ; CALL lbc_lnk( zoty1, ssnd(jps_ocyw)%clgrid, -1. )   
     2291         !   
     2292         !   
     2293         IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components   
     2294         !                                                                        ! Ocean component   
     2295            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->e', ztmp1 )       ! 1st component    
     2296            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->n', ztmp2 )       ! 2nd component    
     2297            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components    
     2298            zoty1(:,:) = ztmp2(:,:)    
     2299            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component   
     2300               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component    
     2301               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component    
     2302               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components    
     2303               zity1(:,:) = ztmp2(:,:)   
     2304            ENDIF   
     2305         ENDIF   
     2306         !   
     2307!         ! spherical coordinates to cartesian -> 2 components to 3 components   
     2308!         IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN   
     2309!            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents   
     2310!            ztmp2(:,:) = zoty1(:,:)   
     2311!            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )   
     2312!            !   
     2313!            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities   
     2314!               ztmp1(:,:) = zitx1(:,:)   
     2315!               ztmp1(:,:) = zity1(:,:)   
     2316!               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )   
     2317!            ENDIF   
     2318!         ENDIF   
     2319         !   
     2320         IF( ssnd(jps_ocxw)%laction )   CALL cpl_snd( jps_ocxw, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid   
     2321         IF( ssnd(jps_ocyw)%laction )   CALL cpl_snd( jps_ocyw, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid   
     2322         !    
     2323      ENDIF   
     2324      !   
     2325      IF( ssnd(jps_ficet)%laction ) THEN   
     2326         CALL cpl_snd( jps_ficet, isec, RESHAPE ( fr_i, (/jpi,jpj,1/) ), info )   
     2327      END IF   
     2328      !                                                      ! ------------------------- !   
     2329      !                                                      !   Water levels to waves   !   
     2330      !                                                      ! ------------------------- !   
     2331      IF( ssnd(jps_wlev)%laction ) THEN   
     2332         IF( ln_apr_dyn ) THEN    
     2333            IF( kt /= nit000 ) THEN    
     2334               ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )    
     2335            ELSE    
     2336               ztmp1(:,:) = sshb(:,:)    
     2337            ENDIF    
     2338         ELSE    
     2339            ztmp1(:,:) = sshn(:,:)    
     2340         ENDIF    
     2341         CALL cpl_snd( jps_wlev  , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )   
     2342      END IF   
    20212343      ! 
    20222344      !  Fields sent by OPA to SAS when doing OPA<->SAS coupling 
  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r7470 r7471  
    8989         &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   & 
    9090         &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   & 
    91          &             nn_lsm    , nn_limflx , nn_components, ln_cpl 
     91         &             ln_tauoc  , ln_stcor  , nn_lsm, nn_limflx , nn_components, ln_cpl 
    9292      INTEGER  ::   ios 
    9393      INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm 
     
    132132         WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl      = ', ln_cpl 
    133133         WRITE(numout,*) '              forced-coupled mixed formulation           ln_mixcpl   = ', ln_mixcpl 
     134         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 
    134139         WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis    = ', lk_oasis 
    135140         WRITE(numout,*) '              components of your executable              nn_components = ', nn_components 
     
    216221       
    217222      IF ( ln_wave ) THEN 
    218       !Activated wave module but neither drag nor stokes drift activated 
    219          IF ( .NOT.(ln_cdgw .OR. ln_sdw) )   THEN 
    220             CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' ) 
    221       !drag coefficient read from wave model definable only with mfs bulk formulae and core  
     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 
    222227         ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN        
    223228             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') 
    224231         ENDIF 
    225232      ELSE 
    226       IF ( ln_cdgw .OR. ln_sdw  )                                         &  
    227          &   CALL ctl_stop('Not Activated Wave Module (ln_wave=F) but     & 
    228          & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ') 
     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)'  ) 
    229239      ENDIF  
    230240      !                          ! Choice of the Surface Boudary Condition (set nsbc) 
     
    382392      END SELECT 
    383393 
     394      IF ( ln_wave .AND. ln_tauoc) THEN                 ! Wave stress subctracted  
     395            utau(:,:) = utau(:,:)*tauoc_wave(:,:)  
     396            vtau(:,:) = vtau(:,:)*tauoc_wave(:,:)  
     397            taum(:,:) = taum(:,:)*tauoc_wave(:,:)  
     398      !  
     399            SELECT CASE( nsbc )  
     400            CASE(  0,1,2,3,5,-1 )  ;  
     401                IF(lwp .AND. kt == nit000 ) WRITE(numout,*) 'WARNING: You are subtracting the wave stress to the ocean. &  
     402                        & If not requested select ln_tauoc=.false'  
     403            END SELECT  
     404      !  
     405      END IF 
    384406      IF( ln_mixcpl )        CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing 
    385407 
  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r7470 r7471  
    44   !! Wave module  
    55   !!====================================================================== 
    6    !! History :  3.3.1  !   2011-09  (Adani M)  Original code: Drag Coefficient  
    7    !!         :  3.4    !   2012-10  (Adani M)                 Stokes Drift  
    8    !!---------------------------------------------------------------------- 
    9    USE iom             ! I/O manager library 
    10    USE in_out_manager  ! I/O manager 
    11    USE lib_mpp         ! distribued memory computing library 
    12    USE fldread        ! read input fields 
    13    USE oce 
    14    USE sbc_oce        ! Surface boundary condition: ocean fields 
    15    USE domvvl 
    16  
    17     
    18    !!---------------------------------------------------------------------- 
    19    !!   sbc_wave       : read drag coefficient from wave model in netcdf files  
    20    !!---------------------------------------------------------------------- 
     6   !! History :  3.3  !   2011-09  (Adani M)  Original code: Drag Coefficient  
     7   !!         :  3.4  !   2012-10  (Adani M)                 Stokes Drift  
     8   !!            3.6  !   2014-09  (Clementi E, Oddo P)New Stokes Drift Computation 
     9   !!---------------------------------------------------------------------- 
     10 
     11   !!---------------------------------------------------------------------- 
     12   !!   sbc_wave      : wave data from wave model in netcdf files  
     13   !!---------------------------------------------------------------------- 
     14   USE oce            !  
     15   USE sbc_oce       ! Surface boundary condition: ocean fields 
     16   USE bdy_oce        ! 
     17   USE domvvl         ! 
     18   ! 
     19   USE iom            ! I/O manager library 
     20   USE in_out_manager ! I/O manager 
     21   USE lib_mpp        ! distribued memory computing library 
     22   USE fldread       ! read input fields 
     23   USE wrk_nemo       ! 
     24   USE phycst         ! physical constants  
    2125 
    2226   IMPLICIT NONE 
    2327   PRIVATE 
    2428 
    25    PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs 
     29   PUBLIC   sbc_stokes, sbc_qiao  ! routines called in sbccpl 
     30   PUBLIC   sbc_wave    ! routine called in sbcmod 
    2631    
    27    INTEGER , PARAMETER ::   jpfld  = 3           ! maximum number of files to read for srokes drift 
    28    INTEGER , PARAMETER ::   jp_usd = 1           ! index of stokes drift  (i-component) (m/s)    at T-point 
    29    INTEGER , PARAMETER ::   jp_vsd = 2           ! index of stokes drift  (j-component) (m/s)    at T-point 
    30    INTEGER , PARAMETER ::   jp_wn  = 3           ! index of wave number                 (1/m)    at T-point 
    31    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
    32    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
    33    REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave  
    34    REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum  
    35    REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d  
     32   ! Variables checking if the wave parameters are coupled (if not, they are read from file) 
     33   LOGICAL, PUBLIC     ::   cpl_hsig=.FALSE. 
     34   LOGICAL, PUBLIC     ::   cpl_phioc=.FALSE. 
     35   LOGICAL, PUBLIC     ::   cpl_sdrftx=.FALSE. 
     36   LOGICAL, PUBLIC     ::   cpl_sdrfty=.FALSE. 
     37   LOGICAL, PUBLIC     ::   cpl_wper=.FALSE. 
     38   LOGICAL, PUBLIC     ::   cpl_wnum=.FALSE. 
     39   LOGICAL, PUBLIC     ::   cpl_wstrf=.FALSE. 
     40   LOGICAL, PUBLIC     ::   cpl_wdrag=.FALSE. 
     41 
     42   INTEGER ::   jpfld                ! number of files to read for stokes drift 
     43   INTEGER ::   jp_usd               ! index of stokes drift  (i-component) (m/s)    at T-point 
     44   INTEGER ::   jp_vsd               ! index of stokes drift  (j-component) (m/s)    at T-point 
     45   INTEGER ::   jp_swh               ! index of significant wave hight      (m)      at T-point 
     46   INTEGER ::   jp_wmp               ! index of mean wave period            (s)      at T-point 
     47 
     48   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
     49   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
     50   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao 
     51   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 
     52   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: cdn_wave  
     53   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: swh,wmp, wnum 
     54   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tauoc_wave 
     55   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tsd2d 
     56   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: zusd2dt, zvsd2dt 
     57   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3d, vsd3d, wsd3d  
     58   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3dt, vsd3dt 
    3659 
    3760   !! * Substitutions 
    3861#  include "domzgr_substitute.h90" 
     62#  include "vectopt_loop_substitute.h90" 
    3963   !!---------------------------------------------------------------------- 
    4064   !! NEMO/OPA 4.0 , NEMO Consortium (2011)  
     
    4468CONTAINS 
    4569 
     70   SUBROUTINE sbc_stokes( ) 
     71      !!--------------------------------------------------------------------- 
     72      !!                     ***  ROUTINE sbc_stokes  *** 
     73      !! 
     74      !! ** Purpose :   compute the 3d Stokes Drift according to Breivik et al., 
     75      !!                2014 (DOI: 10.1175/JPO-D-14-0020.1) 
     76      !! 
     77      !! ** Method  : - Calculate Stokes transport speed  
     78      !!              - Calculate horizontal divergence  
     79      !!              - Integrate the horizontal divergenze from the bottom  
     80      !! ** action   
     81      !!--------------------------------------------------------------------- 
     82      INTEGER                ::   jj,ji,jk  
     83      REAL(wp)                       ::  ztransp, zfac, zsp0, zk, zus, zvs 
     84      REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv   ! 3D workspace 
     85      !!--------------------------------------------------------------------- 
     86      ! 
     87 
     88      CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv ) 
     89      DO jk = 1, jpk 
     90         DO jj = 1, jpj 
     91            DO ji = 1, jpi 
     92               ! On T grid 
     93               ! Stokes transport speed estimated from Hs and Tmean 
     94               ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp)) 
     95               ! Stokes surface speed 
     96               zsp0 = SQRT( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2) 
     97               ! Wavenumber scale 
     98               zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp) 
     99               ! Depth attenuation 
     100               zfac = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk)) 
     101               ! 
     102               usd3dt(ji,jj,jk) = zfac * zusd2dt(ji,jj) * tmask(ji,jj,jk) 
     103               vsd3dt(ji,jj,jk) = zfac * zvsd2dt(ji,jj) * tmask(ji,jj,jk) 
     104            END DO 
     105         END DO 
     106      END DO  
     107      ! Into the U and V Grid 
     108      DO jk = 1, jpkm1 
     109         DO jj = 1, jpjm1 
     110            DO ji = 1, fs_jpim1 
     111               usd3d(ji,jj,jk) = 0.5 *  umask(ji,jj,jk) *   & 
     112                               &  ( usd3dt(ji,jj,jk) + usd3dt(ji+1,jj,jk) ) 
     113               vsd3d(ji,jj,jk) = 0.5 *  vmask(ji,jj,jk) *   & 
     114                               &  ( vsd3dt(ji,jj,jk) + vsd3dt(ji,jj+1,jk) ) 
     115            END DO 
     116         END DO 
     117      END DO 
     118      ! 
     119      CALL lbc_lnk( usd3d(:,:,:), 'U', -1. ) 
     120      CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. ) 
     121      ! 
     122      DO jk = 1, jpkm1               ! Horizontal divergence 
     123         DO jj = 2, jpj 
     124            DO ji = fs_2, jpi 
     125               ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * usd3d(ji  ,jj,jk)     & 
     126                  &                 - e2u(ji-1,jj) * usd3d(ji-1,jj,jk)     & 
     127                  &                 + e1v(ji,jj  ) * vsd3d(ji,jj  ,jk)     & 
     128                  &                 - e1v(ji,jj-1) * vsd3d(ji,jj-1,jk)   ) * r1_e12t(ji,jj) 
     129            END DO 
     130         END DO 
     131      END DO 
     132      ! 
     133      IF( .NOT. AGRIF_Root() ) THEN 
     134         IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   :  ,:) = 0._wp      ! east 
     135         IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   :  ,:) = 0._wp      ! west 
     136         IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :   ,nlcj-1,:) = 0._wp      ! north 
     137         IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2   ,:) = 0._wp      ! south 
     138      ENDIF 
     139      ! 
     140      CALL lbc_lnk( ze3hdiv, 'T', 1. ) 
     141      ! 
     142      DO jk = jpkm1, 1, -1                   ! integrate from the bottom the e3t * hor. divergence 
     143         wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - fse3t_n(:,:,jk) * ze3hdiv(:,:,jk) 
     144      END DO 
     145#if defined key_bdy 
     146      IF( lk_bdy ) THEN 
     147         DO jk = 1, jpkm1 
     148            wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 
     149         END DO 
     150      ENDIF 
     151#endif 
     152      CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv ) 
     153      ! 
     154   END SUBROUTINE sbc_stokes 
     155 
     156   SUBROUTINE sbc_qiao 
     157      !!--------------------------------------------------------------------- 
     158      !!                     ***  ROUTINE sbc_qiao  *** 
     159      !! 
     160      !! ** Purpose :   Qiao formulation for wave enhanced turbulence 
     161      !!                2010 (DOI: 10.1007/s10236-010-0326)  
     162      !! 
     163      !! ** Method  : -  
     164      !! ** action   
     165      !!--------------------------------------------------------------------- 
     166      INTEGER :: jj, ji 
     167 
     168      ! Calculate the module of the stokes drift on T grid 
     169      !------------------------------------------------- 
     170      DO jj = 1, jpj 
     171         DO ji = 1, jpi 
     172            tsd2d(ji,jj) = SQRT( zusd2dt(ji,jj) * zusd2dt(ji,jj) + zvsd2dt(ji,jj) * zvsd2dt(ji,jj) ) 
     173         END DO 
     174      END DO 
     175      ! 
     176   END SUBROUTINE sbc_qiao 
     177 
    46178   SUBROUTINE sbc_wave( kt ) 
    47179      !!--------------------------------------------------------------------- 
    48       !!                     ***  ROUTINE sbc_apr  *** 
    49       !! 
    50       !! ** Purpose :   read drag coefficient from wave model  in netcdf files. 
     180      !!                     ***  ROUTINE sbc_wave  *** 
     181      !! 
     182      !! ** Purpose :   read wave parameters from wave model  in netcdf files. 
    51183      !! 
    52184      !! ** Method  : - Read namelist namsbc_wave 
    53185      !!              - Read Cd_n10 fields in netcdf files  
    54186      !!              - Read stokes drift 2d in netcdf files  
    55       !!              - Read wave number      in netcdf files  
    56       !!              - Compute 3d stokes drift using monochromatic 
    57       !! ** action  :    
    58       !!                
    59       !!--------------------------------------------------------------------- 
    60       USE oce,  ONLY : un,vn,hdivn,rotn 
    61       USE divcur 
    62       USE wrk_nemo 
    63 #if defined key_bdy 
    64       USE bdy_oce, ONLY : bdytmask 
    65 #endif 
    66       INTEGER, INTENT( in  ) ::  kt       ! ocean time step 
    67       INTEGER                ::  ierror   ! return error code 
    68       INTEGER                ::  ifpr, jj,ji,jk  
    69       INTEGER                ::   ios     ! Local integer output status for namelist read 
    70       REAL(wp),DIMENSION(:,:,:),POINTER             ::  udummy,vdummy,hdivdummy,rotdummy 
    71       REAL                                          ::  z2dt,z1_2dt 
    72       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read 
     187      !!              - Read wave number in netcdf files  
     188      !!              - Compute 3d stokes drift using Breivik et al.,2014 
     189      !!                formulation 
     190      !! ** action   
     191      !!--------------------------------------------------------------------- 
     192      USE zdf_oce,  ONLY : ln_zdfqiao 
     193 
     194      IMPLICIT NONE 
     195 
     196      INTEGER, INTENT( in  ) ::   kt       ! ocean time step 
     197      ! 
     198      INTEGER                ::   ierror   ! return error code 
     199      INTEGER                ::   ifpr 
     200      INTEGER                ::   ios      ! Local integer output status for namelist read 
     201      ! 
    73202      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files 
    74       TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_wn   ! informations about the fields to be read 
    75       !!--------------------------------------------------------------------- 
    76       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn 
    77       !!--------------------------------------------------------------------- 
    78  
    79       !!---------------------------------------------------------------------- 
    80       ! 
     203      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i     ! array of namelist informations on the fields to read 
     204      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  & 
     205                             &   sn_swh, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read 
     206      !! 
     207      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_swh, sn_wmp, sn_wnum, sn_tauoc 
     208      !!--------------------------------------------------------------------- 
    81209      ! 
    82210      !                                         ! -------------------- ! 
     
    92220         IF(lwm) WRITE ( numond, namsbc_wave ) 
    93221         ! 
    94  
    95          IF ( ln_cdgw ) THEN 
    96             ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    97             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    98             ! 
    99                                    ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
    100             IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
    101             CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
     222         IF( ln_cdgw ) THEN 
     223            IF( .NOT. cpl_wdrag ) THEN 
     224               ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
     225               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     226               ! 
     227                                      ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
     228               IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
     229               CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
     230            ENDIF 
    102231            ALLOCATE( cdn_wave(jpi,jpj) ) 
    103232            cdn_wave(:,:) = 0.0 
    104         ENDIF 
    105          IF ( ln_sdw ) THEN 
    106             slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn 
    107             ALLOCATE( sf_sd(3), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    108             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    109             ! 
    110             DO ifpr= 1, jpfld 
    111                ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
    112                IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
    113             END DO 
    114             CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
    115             ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj),uwavenum(jpi,jpj),vwavenum(jpi,jpj) ) 
     233         ENDIF 
     234 
     235         IF( ln_tauoc ) THEN 
     236            IF( .NOT. cpl_wstrf ) THEN 
     237               ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc 
     238               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     239               ! 
     240                                       ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   ) 
     241               IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 
     242               CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
     243            ENDIF 
     244            ALLOCATE( tauoc_wave(jpi,jpj) ) 
     245         ENDIF 
     246 
     247         IF( ln_sdw ) THEN 
     248            ! Find out how many fields have to be read from file if not coupled 
     249            jpfld=0 
     250            jp_usd=0; jp_vsd=0; jp_swh=0; jp_wmp=0 
     251            IF( .NOT. cpl_sdrftx ) THEN 
     252               jpfld=jpfld+1 
     253               jp_usd=jpfld 
     254            ENDIF 
     255            IF( .NOT. cpl_sdrfty ) THEN 
     256               jpfld=jpfld+1 
     257               jp_vsd=jpfld 
     258            ENDIF 
     259            IF( .NOT. cpl_hsig ) THEN 
     260               jpfld=jpfld+1 
     261               jp_swh=jpfld 
     262            ENDIF 
     263            IF( .NOT. cpl_wper ) THEN 
     264               jpfld=jpfld+1 
     265               jp_wmp=jpfld 
     266            ENDIF 
     267 
     268            ! Read from file only the non-coupled fields  
     269            IF( jpfld > 0 ) THEN 
     270               ALLOCATE( slf_i(jpfld) ) 
     271               IF( jp_usd > 0 ) slf_i(jp_usd) = sn_usd 
     272               IF( jp_vsd > 0 ) slf_i(jp_vsd) = sn_vsd 
     273               IF( jp_swh > 0 ) slf_i(jp_swh) = sn_swh 
     274               IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 
     275               ALLOCATE( sf_sd(jpfld), STAT=ierror )           !* allocate and fill sf_sd with stokes drift 
     276               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     277               ! 
     278               DO ifpr= 1, jpfld 
     279                  ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
     280                  IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
     281               END DO 
     282 
     283               CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
     284            ENDIF 
    116285            ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) ) 
    117             usd2d(:,:) = 0.0 ;  vsd2d(:,:) = 0.0 ; uwavenum(:,:) = 0.0 ; vwavenum(:,:) = 0.0 
    118             usd3d(:,:,:) = 0.0 ;vsd3d(:,:,:) = 0.0 ; wsd3d(:,:,:) = 0.0 
    119          ENDIF 
    120       ENDIF 
     286            ALLOCATE( usd3dt(jpi,jpj,jpk),vsd3dt(jpi,jpj,jpk) ) 
     287            ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) ) 
     288            ALLOCATE( zusd2dt(jpi,jpj), zvsd2dt(jpi,jpj) ) 
     289            usd3d(:,:,:) = 0._wp 
     290            vsd3d(:,:,:) = 0._wp 
     291            wsd3d(:,:,:) = 0._wp 
     292            IF( ln_zdfqiao ) THEN     !==  Vertical mixing enhancement using Qiao,2010  ==! 
     293               IF( .NOT. cpl_wnum ) THEN 
     294                  ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum 
     295                  IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable toallocate sf_wave structure' ) 
     296                                         ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
     297                  IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
     298                  CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
     299               ENDIF 
     300               ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) ) 
     301            ENDIF 
     302         ENDIF 
     303      ENDIF 
     304      ! 
     305      IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN              !==  Neutral drag coefficient  ==! 
     306         CALL fld_read( kt, nn_fsbc, sf_cd )      ! read from external forcing 
     307         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
     308      ENDIF 
     309 
     310      IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN             !==  Wave induced stress  ==! 
     311         CALL fld_read( kt, nn_fsbc, sf_tauoc )      !* read wave norm stress from external forcing 
     312         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 
     313      ENDIF 
     314 
     315      IF( ln_sdw )  THEN                         !==  Computation of the 3d Stokes Drift  ==!  
    121316         ! 
     317         ! Read from file only if the field is not coupled 
     318         IF( jpfld > 0 ) THEN 
     319            CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing 
     320            IF( jp_swh > 0 ) swh(:,:)     = sf_sd(jp_swh)%fnow(:,:,1)   ! significant wave height 
     321            IF( jp_wmp > 0 ) wmp(:,:)     = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
     322            IF( jp_usd > 0 ) zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
     323            IF( jp_vsd > 0 ) zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
     324         ENDIF 
    122325         ! 
    123       IF ( ln_cdgw ) THEN 
    124          CALL fld_read( kt, nn_fsbc, sf_cd )      !* read drag coefficient from external forcing 
    125          cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
    126       ENDIF 
    127       IF ( ln_sdw )  THEN 
    128           CALL fld_read( kt, nn_fsbc, sf_sd )      !* read drag coefficient from external forcing 
    129  
    130          ! Interpolate wavenumber, stokes drift into the grid_V and grid_V 
    131          !------------------------------------------------- 
    132  
    133          DO jj = 1, jpjm1 
    134             DO ji = 1, jpim1 
    135                uwavenum(ji,jj)=0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    136                &                                + sf_sd(3)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) ) 
    137  
    138                vwavenum(ji,jj)=0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    139                &                                + sf_sd(3)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) ) 
    140  
    141                usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(1)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    142                &                                + sf_sd(1)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) ) 
    143  
    144                vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(2)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    145                &                                + sf_sd(2)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) ) 
    146             END DO 
    147          END DO 
    148  
    149           !Computation of the 3d Stokes Drift 
    150           DO jk = 1, jpk 
    151              DO jj = 1, jpj-1 
    152                 DO ji = 1, jpi-1 
    153                    usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk)))) 
    154                    vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk)))) 
    155                 END DO 
    156              END DO 
    157              usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) ) 
    158              vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,jpj,jk)) ) 
    159           END DO 
    160  
    161           CALL wrk_alloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 
    162            
    163           udummy(:,:,:)=un(:,:,:) 
    164           vdummy(:,:,:)=vn(:,:,:) 
    165           hdivdummy(:,:,:)=hdivn(:,:,:) 
    166           rotdummy(:,:,:)=rotn(:,:,:) 
    167           un(:,:,:)=usd3d(:,:,:) 
    168           vn(:,:,:)=vsd3d(:,:,:) 
    169           CALL div_cur(kt) 
    170       !                                           !------------------------------! 
    171       !                                           !     Now Vertical Velocity    ! 
    172       !                                           !------------------------------! 
    173           z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
    174  
    175           z1_2dt = 1.e0 / z2dt 
    176           DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence 
    177              ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 
    178              wsd3d(:,:,jk) = wsd3d(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        & 
    179                 &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    & 
    180                 &                         * tmask(:,:,jk) * z1_2dt 
    181 #if defined key_bdy 
    182              wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 
    183 #endif 
    184           END DO 
    185           hdivn(:,:,:)=hdivdummy(:,:,:) 
    186           rotn(:,:,:)=rotdummy(:,:,:) 
    187           vn(:,:,:)=vdummy(:,:,:) 
    188           un(:,:,:)=udummy(:,:,:) 
    189           CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 
    190       ENDIF 
     326         ! Read also wave number if needed, so that it is available in coupling routines 
     327         IF( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 
     328            CALL fld_read( kt, nn_fsbc, sf_wn )      !* read wave parameters from external forcing 
     329            wnum(:,:) = sf_wn(1)%fnow(:,:,1) 
     330         ENDIF 
     331            
     332         !==  Computation of the 3d Stokes Drift according to Breivik et al.,2014 
     333         !(DOI: 10.1175/JPO-D-14-0020.1)==!  
     334         ! 
     335         ! Calculate only if no necessary fields are coupled, if not calculate later after coupling 
     336         IF( jpfld == 4 ) THEN 
     337            CALL sbc_stokes() 
     338            IF( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 
     339               CALL sbc_qiao() 
     340            ENDIF 
     341         ENDIF 
     342      ENDIF 
     343      ! 
    191344   END SUBROUTINE sbc_wave 
    192345       
  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

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

    r7470 r7471  
    3535   INTEGER , PUBLIC ::   nn_npc      !: non penetrative convective scheme call  frequency 
    3636   INTEGER , PUBLIC ::   nn_npcp     !: non penetrative convective scheme print frequency 
     37   LOGICAL , PUBLIC ::   ln_zdfqiao  !: Enhanced wave vertical mixing Qiao(2010) formulation flag 
    3738 
    3839 
  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90

    r7470 r7471  
    5353      INTEGER ::   ios 
    5454      !! 
    55       NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,   & 
    56          &              ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp 
     55      NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,  &  
     56         &        ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp,       &  
     57         &        ln_zdfqiao  
    5758      !!---------------------------------------------------------------------- 
    5859 
     
    8384         WRITE(numout,*) '      npc call  frequency                 nn_npc    = ', nn_npc 
    8485         WRITE(numout,*) '      npc print frequency                 nn_npcp   = ', nn_npcp 
     86         WRITE(numout,*) '      Qiao formulation flag               ln_zdfqiao=', ln_zdfqiao 
    8587      ENDIF 
    8688 
  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step.F90

    r7470 r7471  
    2525   !!            3.4  !  2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal 
    2626   !!                 !  2012-07  (J. Simeon, G. Madec, C. Ethe) Online coarsening of outputs 
     27   !!            3.6  !  2014-10  (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves 
    2728   !!            3.7  !  2014-04  (F. Roquet, G. Madec) New equations of state 
    2829   !!---------------------------------------------------------------------- 
     
    7273      !!              -8- Outputs and diagnostics 
    7374      !!---------------------------------------------------------------------- 
    74       INTEGER ::   jk      ! dummy loop indice 
     75      INTEGER ::   ji,jj,jk ! dummy loop indice 
    7576      INTEGER ::   indic    ! error indicator if < 0 
    7677      INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt) 
     
    132133      IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz 
    133134      IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz 
     135      IF( ln_zdfqiao )   CALL zdf_qiao( kstp )           ! Qiao vertical mixing   
     136      !  
    134137      IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz 
    135138      IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value) 
     
    209212                                  CALL dyn_adv      ( kstp )   ! advection (vector or flux form) 
    210213                                  CALL dyn_vor      ( kstp )   ! vorticity term including Coriolis 
     214          IF( ln_wave .AND. ln_sdw .AND. ln_stcor )          &  
     215                        &         CALL dyn_stcor     ( kstp )  ! Stokes-Coriolis forcing 
    211216                                  CALL dyn_ldf      ( kstp )   ! lateral mixing 
    212217          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! add Neptune velocities (simplified) 
  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r7470 r7471  
    2828   USE sbctide          ! Tide initialisation 
    2929   USE sbcapr           ! surface boundary condition: ssh_ib required by bdydta  
     30   USE sbcwave          ! Wave intialisation 
    3031 
    3132   USE traqsr           ! solar radiation penetration      (tra_qsr routine) 
     
    5152   USE dynspg_oce       ! surface pressure gradient        (dyn_spg routine) 
    5253   USE dynspg           ! surface pressure gradient        (dyn_spg routine) 
     54   USE dynstcor         ! simp. form of Stokes-Coriolis 
    5355   USE dynnept          ! simp. form of Neptune effect(dyn_nept_cor routine) 
    5456 
     
    8486   USE zdfric           ! Richardson vertical mixing       (zdf_ric routine) 
    8587   USE zdfmxl           ! Mixed-layer depth                (zdf_mxl routine) 
     88   USE zdfqiao          !Qiao module wave induced mixing   (zdf_qiao routine) 
    8689 
    8790   USE zpshde           ! partial step: hor. derivative     (zps_hde routine) 
  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/SETTE/sette.sh

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