New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 14066 for NEMO/branches/2020/dev_r13787_doc_latex_recovery/src/OCE/SBC/sbcmod.F90 – NEMO

Ignore:
Timestamp:
2020-12-03T20:14:08+01:00 (4 years ago)
Author:
nicolasmartin
Message:

#2414 Sync merge with trunk

Location:
NEMO/branches/2020/dev_r13787_doc_latex_recovery
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13787_doc_latex_recovery

    • Property svn:externals
      •  

        old new  
        88 
        99# SETTE 
        10 ^/utils/CI/sette@13559        sette 
         10^/utils/CI/sette_wave@13990         sette 
  • NEMO/branches/2020/dev_r13787_doc_latex_recovery/src/OCE/SBC/sbcmod.F90

    r13722 r14066  
    1616   !!            4.0  ! 2016-06  (L. Brodeau) new general bulk formulation 
    1717   !!            4.0  ! 2019-03  (F. Lemarié & G. Samson)  add ABL compatibility (ln_abl=TRUE) 
     18   !!            4.2  ! 2020-12  (G. Madec, E. Clementi) modified wave forcing and coupling   
    1819   !!---------------------------------------------------------------------- 
    1920 
     
    5455   USE usrdef_sbc     ! user defined: surface boundary condition 
    5556   USE closea         ! closed sea 
     57   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    5658   ! 
    5759   USE prtctl         ! Print control                    (prt_ctl routine) 
     
    7072 
    7173   INTEGER ::   nsbc   ! type of surface boundary condition (deduced from namsbc informations) 
    72  
     74   !! * Substitutions 
     75#  include "do_loop_substitute.h90" 
    7376   !!---------------------------------------------------------------------- 
    7477   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    99102         &             nn_ice   , ln_ice_embd,                                       & 
    100103         &             ln_traqsr, ln_dm2dc ,                                         & 
    101          &             ln_rnf   , nn_fwb   , ln_ssr   , ln_apr_dyn,                  & 
    102          &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauwoc , ln_stcor  ,      & 
    103          &             ln_tauw  , nn_lsm, nn_sdrift 
     104         &             ln_rnf   , nn_fwb     , ln_ssr   , ln_apr_dyn,                & 
     105         &             ln_wave  , nn_lsm 
    104106      !!---------------------------------------------------------------------- 
    105107      ! 
     
    133135         WRITE(numout,*) '         bulk         formulation                   ln_blk        = ', ln_blk 
    134136         WRITE(numout,*) '         ABL          formulation                   ln_abl        = ', ln_abl 
     137         WRITE(numout,*) '         Surface wave (forced or coupled)           ln_wave       = ', ln_wave 
    135138         WRITE(numout,*) '      Type of coupling (Ocean/Ice/Atmosphere) : ' 
    136139         WRITE(numout,*) '         ocean-atmosphere coupled formulation       ln_cpl        = ', ln_cpl 
     
    150153         WRITE(numout,*) '         runoff / runoff mouths                     ln_rnf        = ', ln_rnf 
    151154         WRITE(numout,*) '         nb of iterations if land-sea-mask applied  nn_lsm        = ', nn_lsm 
    152          WRITE(numout,*) '         surface wave                               ln_wave       = ', ln_wave 
    153          WRITE(numout,*) '               Stokes drift corr. to vert. velocity ln_sdw        = ', ln_sdw 
    154          WRITE(numout,*) '                  vertical parametrization          nn_sdrift     = ', nn_sdrift 
    155          WRITE(numout,*) '               wave modified ocean stress           ln_tauwoc     = ', ln_tauwoc 
    156          WRITE(numout,*) '               wave modified ocean stress component ln_tauw       = ', ln_tauw 
    157          WRITE(numout,*) '               Stokes coriolis term                 ln_stcor      = ', ln_stcor 
    158          WRITE(numout,*) '               neutral drag coefficient (CORE,NCAR) ln_cdgw       = ', ln_cdgw 
    159       ENDIF 
    160       ! 
    161       IF( .NOT.ln_wave ) THEN 
    162          ln_sdw = .false. ; ln_cdgw = .false. ; ln_tauwoc = .false. ; ln_tauw = .false. ; ln_stcor = .false. 
    163       ENDIF  
    164       IF( ln_sdw ) THEN 
    165          IF( .NOT.(nn_sdrift==jp_breivik_2014 .OR. nn_sdrift==jp_li_2017 .OR. nn_sdrift==jp_peakfr) ) & 
    166             CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' ) 
    167       ENDIF 
    168       ll_st_bv2014  = ( nn_sdrift==jp_breivik_2014 ) 
    169       ll_st_li2017  = ( nn_sdrift==jp_li_2017 ) 
    170       ll_st_bv_li   = ( ll_st_bv2014 .OR. ll_st_li2017 ) 
    171       ll_st_peakfr  = ( nn_sdrift==jp_peakfr ) 
    172       IF( ln_tauwoc .AND. ln_tauw ) & 
    173          CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', & 
    174                                   '(ln_tauwoc=.true. and ln_tauw=.true.)' ) 
    175       IF( ln_tauwoc ) & 
    176          CALL ctl_warn( 'You are subtracting the wave stress to the ocean (ln_tauwoc=.true.)' ) 
    177       IF( ln_tauw ) & 
    178          CALL ctl_warn( 'The wave modified ocean stress components are used (ln_tauw=.true.) ', & 
    179                               'This will override any other specification of the ocean stress' ) 
     155      ENDIF 
    180156      ! 
    181157      IF( .NOT.ln_usr ) THEN     ! the model calendar needs some specificities (except in user defined case) 
     
    357333      IF( nn_ice == 3 )   CALL cice_sbc_init( nsbc, Kbb, Kmm )   ! CICE initialization 
    358334      ! 
    359       IF( ln_wave     )   CALL sbc_wave_init                     ! surface wave initialisation 
    360       ! 
    361       IF( lwxios ) THEN 
    362          CALL iom_set_rstw_var_active('utau_b') 
    363          CALL iom_set_rstw_var_active('vtau_b') 
    364          CALL iom_set_rstw_var_active('qns_b') 
    365          ! The 3D heat content due to qsr forcing is treated in traqsr 
    366          ! CALL iom_set_rstw_var_active('qsr_b') 
    367          CALL iom_set_rstw_var_active('emp_b') 
    368          CALL iom_set_rstw_var_active('sfx_b') 
    369       ENDIF 
    370  
     335      IF( ln_wave     ) THEN 
     336                          CALL sbc_wave_init                     ! surface wave initialisation 
     337      ELSE 
     338                          IF(lwp) WRITE(numout,*) 
     339                          IF(lwp) WRITE(numout,*) '   No surface waves : all wave related logical set to false' 
     340                          ln_sdw       = .false. 
     341                          ln_stcor     = .false. 
     342                          ln_cdgw      = .false. 
     343                          ln_tauoc     = .false. 
     344                          ln_wave_test = .false. 
     345                          ln_charn     = .false. 
     346                          ln_taw       = .false. 
     347                          ln_phioc     = .false. 
     348                          ln_bern_srfc = .false. 
     349                          ln_breivikFV_2016 = .false. 
     350                          ln_vortex_force = .false. 
     351                          ln_stshear  = .false. 
     352      ENDIF 
     353      ! 
    371354   END SUBROUTINE sbc_init 
    372355 
     
    390373      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    391374      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices 
     375      INTEGER  ::   jj, ji          ! dummy loop argument 
    392376      ! 
    393377      LOGICAL ::   ll_sas, ll_opa   ! local logical 
     
    422406      ! 
    423407      IF( .NOT.ll_sas )   CALL sbc_ssm ( kt, Kbb, Kmm )  ! mean ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m) 
    424       IF( ln_wave     )   CALL sbc_wave( kt, Kmm )       ! surface waves 
    425  
    426408      ! 
    427409      !                                            !==  sbc formulation  ==! 
    428410      !                                                    
     411      ! 
    429412      SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition 
    430413      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx) 
     
    433416      CASE( jp_blk     ) 
    434417         IF( ll_sas    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )   ! OPA-SAS coupling: SAS receiving fields from OPA 
     418!!!!!!!!!!! ATTENTION:ln_wave is not only used for oasis coupling !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     419         IF( ln_wave )   THEN 
     420             IF ( lk_oasis )  CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm )   ! OPA-wave coupling 
     421             CALL sbc_wave ( kt, Kmm ) 
     422         ENDIF 
    435423                               CALL sbc_blk       ( kt )                    ! bulk formulation for the ocean 
    436424                               ! 
     
    446434      IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! forced-coupled mixed formulation after forcing 
    447435      ! 
    448       IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )              ! Wind stress provided by waves  
     436      IF( ln_wave .AND. ln_tauoc )  THEN            ! Wave stress reduction 
     437         DO_2D( 0, 0, 0, 0) 
     438            utau(ji,jj) = utau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji-1,jj) ) * 0.5_wp 
     439            vtau(ji,jj) = vtau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji,jj-1) ) * 0.5_wp 
     440         END_2D 
     441         ! 
     442         CALL lbc_lnk( 'sbcwave', utau, 'U', -1. ) 
     443         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. ) 
     444         ! 
     445         taum(:,:) = taum(:,:)*tauoc_wave(:,:) 
     446         ! 
     447         IF( kt == nit000 )   CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.',   & 
     448            &                                'If not requested select ln_tauoc=.false.' ) 
     449         ! 
     450      ELSEIF( ln_wave .AND. ln_taw ) THEN                  ! Wave stress reduction 
     451         utau(:,:) = utau(:,:) - tawx(:,:) + twox(:,:) 
     452         vtau(:,:) = vtau(:,:) - tawy(:,:) + twoy(:,:) 
     453         CALL lbc_lnk( 'sbcwave', utau, 'U', -1. ) 
     454         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. ) 
     455         ! 
     456         DO_2D( 0, 0, 0, 0) 
     457             taum(ji,jj) = sqrt((.5*(utau(ji-1,jj)+utau(ji,jj)))**2 + (.5*(vtau(ji,jj-1)+vtau(ji,jj)))**2) 
     458         END_2D 
     459         ! 
     460         IF( kt == nit000 )   CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.',   & 
     461            &                                'If not requested select ln_taw=.false.' ) 
     462         ! 
     463      ENDIF 
     464      CALL lbc_lnk( 'sbcmod', taum(:,:), 'T', 1. ) 
    449465      ! 
    450466      !                                            !==  Misc. Options  ==! 
     
    459475 
    460476      IF( ln_icebergs    )   THEN 
    461                                      CALL icb_stp( kt )           ! compute icebergs 
     477                                     CALL icb_stp( kt, Kmm )           ! compute icebergs 
    462478         ! Icebergs do not melt over the haloes.  
    463479         ! So emp values over the haloes are no more consistent with the inner domain values.  
     
    507523      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    ! 
    508524         !                                             ! ---------------------------------------- ! 
    509          IF( ln_rstart .AND.    &                               !* Restart: read in restart file 
    510             & iom_varid( numror, 'utau_b', ldstop = .FALSE. ) > 0 ) THEN 
    511             IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file' 
    512             CALL iom_get( numror, jpdom_auto, 'utau_b', utau_b, ldxios = lrxios, cd_type = 'U', psgn = -1._wp )   ! before i-stress  (U-point) 
    513             CALL iom_get( numror, jpdom_auto, 'vtau_b', vtau_b, ldxios = lrxios, cd_type = 'V', psgn = -1._wp )   ! before j-stress  (V-point) 
    514             CALL iom_get( numror, jpdom_auto,  'qns_b',  qns_b, ldxios = lrxios )   ! before non solar heat flux (T-point) 
    515             ! The 3D heat content due to qsr forcing is treated in traqsr 
    516             ! CALL iom_get( numror, jpdom_auto, 'qsr_b' , qsr_b, ldxios = lrxios  ) ! before     solar heat flux (T-point) 
    517             CALL iom_get( numror, jpdom_auto, 'emp_b', emp_b, ldxios = lrxios  )    ! before     freshwater flux (T-point) 
     525         IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN            !* Restart: read in restart file 
     526            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields read in the restart file' 
     527            CALL iom_get( numror, jpdom_auto, 'utau_b', utau_b )   ! i-stress 
     528            CALL iom_get( numror, jpdom_auto, 'vtau_b', vtau_b )   ! j-stress 
     529            CALL iom_get( numror, jpdom_auto,  'qns_b',  qns_b )   ! non solar heat flux 
     530            CALL iom_get( numror, jpdom_auto,  'emp_b',  emp_b )   ! freshwater flux 
     531            ! NB: The 3D heat content due to qsr forcing (qsr_hc_b) is treated in traqsr 
    518532            ! To ensure restart capability with 3.3x/3.4 restart files    !! to be removed in v3.6 
    519533            IF( iom_varid( numror, 'sfx_b', ldstop = .FALSE. ) > 0 ) THEN 
    520                CALL iom_get( numror, jpdom_auto, 'sfx_b', sfx_b, ldxios = lrxios )  ! before salt flux (T-point) 
     534               CALL iom_get( numror, jpdom_auto, 'sfx_b', sfx_b )  ! before salt flux (T-point) 
    521535            ELSE 
    522536               sfx_b (:,:) = sfx(:,:) 
     
    538552            &                    'at it= ', kt,' date= ', ndastp 
    539553         IF(lwp) WRITE(numout,*) '~~~~' 
    540          IF( lwxios ) CALL iom_swap(      cwxios_context          ) 
    541          CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau, ldxios = lwxios ) 
    542          CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau, ldxios = lwxios ) 
    543          CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns, ldxios = lwxios  ) 
     554         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau ) 
     555         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau ) 
     556         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  ) 
    544557         ! The 3D heat content due to qsr forcing is treated in traqsr 
    545558         ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  ) 
    546          CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp, ldxios = lwxios  ) 
    547          CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx, ldxios = lwxios  ) 
    548          IF( lwxios ) CALL iom_swap(      cxios_context          ) 
     559         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  ) 
     560         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx  ) 
    549561      ENDIF 
    550562      !                                                ! ---------------------------------------- ! 
     
    552564      !                                                ! ---------------------------------------- ! 
    553565      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
    554          CALL iom_put( "empmr"  , emp    - rnf )                ! upward water flux 
    555          CALL iom_put( "empbmr" , emp_b  - rnf )                ! before upward water flux ( needed to recalculate the time evolution of ssh in offline ) 
    556          CALL iom_put( "saltflx", sfx  )                        ! downward salt flux (includes virtual salt flux beneath ice in linear free surface case) 
    557          CALL iom_put( "fmmflx", fmmflx  )                      ! Freezing-melting water flux 
    558          CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux 
    559          CALL iom_put( "qns"   , qns        )                   ! solar heat flux 
    560          CALL iom_put( "qsr"   ,       qsr  )                   ! solar heat flux 
     566         CALL iom_put( "empmr"  , emp   - rnf )                ! upward water flux 
     567         CALL iom_put( "empbmr" , emp_b - rnf )                ! before upward water flux ( needed to recalculate the time evolution of ssh in offline ) 
     568         CALL iom_put( "saltflx", sfx         )                ! downward salt flux (includes virtual salt flux beneath ice in linear free surface case) 
     569         CALL iom_put( "fmmflx" , fmmflx      )                ! Freezing-melting water flux 
     570         CALL iom_put( "qt"     , qns + qsr   )                ! total heat flux 
     571         CALL iom_put( "qns"    , qns         )                ! solar heat flux 
     572         CALL iom_put( "qsr"    ,       qsr   )                ! solar heat flux 
    561573         IF( nn_ice > 0 .OR. ll_opa )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction 
    562          CALL iom_put( "taum"  , taum       )                   ! wind stress module 
    563          CALL iom_put( "wspd"  , wndm       )                   ! wind speed  module over free ocean or leads in presence of sea-ice 
    564          CALL iom_put( "qrp", qrp )                             ! heat flux damping 
    565          CALL iom_put( "erp", erp )                             ! freshwater flux damping 
     574         CALL iom_put( "taum"   , taum        )                ! wind stress module 
     575         CALL iom_put( "wspd"   , wndm        )                ! wind speed  module over free ocean or leads in presence of sea-ice 
     576         CALL iom_put( "qrp"    , qrp         )                ! heat flux damping 
     577         CALL iom_put( "erp"    , erp         )                ! freshwater flux damping 
    566578      ENDIF 
    567579      ! 
    568580      IF(sn_cfctl%l_prtctl) THEN     ! print mean trends (used for debugging) 
    569          CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i     - : ', mask1=tmask ) 
    570          CALL prt_ctl(tab2d_1=(emp-rnf)        , clinfo1=' emp-rnf  - : ', mask1=tmask ) 
    571          CALL prt_ctl(tab2d_1=(sfx-rnf)        , clinfo1=' sfx-rnf  - : ', mask1=tmask ) 
    572          CALL prt_ctl(tab2d_1=qns              , clinfo1=' qns      - : ', mask1=tmask ) 
    573          CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr      - : ', mask1=tmask ) 
    574          CALL prt_ctl(tab3d_1=tmask            , clinfo1=' tmask    - : ', mask1=tmask, kdim=jpk ) 
     581         CALL prt_ctl(tab2d_1=fr_i                , clinfo1=' fr_i     - : ', mask1=tmask ) 
     582         CALL prt_ctl(tab2d_1=(emp-rnf)           , clinfo1=' emp-rnf  - : ', mask1=tmask ) 
     583         CALL prt_ctl(tab2d_1=(sfx-rnf)           , clinfo1=' sfx-rnf  - : ', mask1=tmask ) 
     584         CALL prt_ctl(tab2d_1=qns                 , clinfo1=' qns      - : ', mask1=tmask ) 
     585         CALL prt_ctl(tab2d_1=qsr                 , clinfo1=' qsr      - : ', mask1=tmask ) 
     586         CALL prt_ctl(tab3d_1=tmask               , clinfo1=' tmask    - : ', mask1=tmask, kdim=jpk ) 
    575587         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' sst      - : ', mask1=tmask, kdim=1   ) 
    576588         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sss      - : ', mask1=tmask, kdim=1   ) 
Note: See TracChangeset for help on using the changeset viewer.