source: NEMO/branches/UKMO/dev_10448_WAD_SBC_BUGFIX/src/OCE/SBC/sbcmod.F90 @ 10456

Last change on this file since 10456 was 10456, checked in by deazer, 22 months ago

Added option to taper sbc fluxes near very shallow water when using WAD
Corrected some IO bugs in dia25h, diatmb for WAD case.
User has control of the tapering. At what depth to start it, and at what fraction to start
the tanh tapering. At the WAD limit SBC is turned off completely.
Dry cells do not have any communication with the atmosphere
To DO: Documentation update.
Although not all sette tests are passed (AGRIF etc.)
it does no worse than the trunk at the revision the branch is made

  • Property svn:keywords set to Id
File size: 34.2 KB
Line 
1MODULE sbcmod
2   !!======================================================================
3   !!                       ***  MODULE  sbcmod  ***
4   !! Surface module :  provide to the ocean its surface boundary condition
5   !!======================================================================
6   !! History :  3.0  ! 2006-07  (G. Madec)  Original code
7   !!            3.1  ! 2008-08  (S. Masson, A. Caubel, E. Maisonnave, G. Madec) coupled interface
8   !!            3.3  ! 2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
9   !!            3.3  ! 2010-10  (S. Masson)  add diurnal cycle
10   !!            3.3  ! 2010-09  (D. Storkey) add ice boundary conditions (BDY)
11   !!             -   ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
12   !!             -   ! 2010-10  (J. Chanut, C. Bricaud, G. Madec)  add the surface pressure forcing
13   !!            3.4  ! 2011-11  (C. Harris) CICE added as an option
14   !!            3.5  ! 2012-11  (A. Coward, G. Madec) Rethink of heat, mass and salt surface fluxes
15   !!            3.6  ! 2014-11  (P. Mathiot, C. Harris) add ice shelves melting
16   !!            4.0  ! 2016-06  (L. Brodeau) new general bulk formulation
17   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   sbc_init      : read namsbc namelist
21   !!   sbc           : surface ocean momentum, heat and freshwater boundary conditions
22   !!   sbc_final     : Finalize CICE ice model (if used)
23   !!----------------------------------------------------------------------
24   USE oce            ! ocean dynamics and tracers
25   USE dom_oce        ! ocean space and time domain
26   USE phycst         ! physical constants
27   USE sbc_oce        ! Surface boundary condition: ocean fields
28   USE trc_oce        ! shared ocean-passive tracers variables
29   USE sbc_ice        ! Surface boundary condition: ice fields
30   USE sbcdcy         ! surface boundary condition: diurnal cycle
31   USE sbcssm         ! surface boundary condition: sea-surface mean variables
32   USE sbcflx         ! surface boundary condition: flux formulation
33   USE sbcblk         ! surface boundary condition: bulk formulation
34   USE sbcice_if      ! surface boundary condition: ice-if sea-ice model
35#if defined key_si3
36   USE icestp         ! surface boundary condition: SI3 sea-ice model
37#endif
38   USE sbcice_cice    ! surface boundary condition: CICE sea-ice model
39   USE sbcisf         ! surface boundary condition: ice-shelf
40   USE sbccpl         ! surface boundary condition: coupled formulation
41   USE cpl_oasis3     ! OASIS routines for coupling
42   USE sbcssr         ! surface boundary condition: sea surface restoring
43   USE sbcrnf         ! surface boundary condition: runoffs
44   USE sbcapr         ! surface boundary condition: atmo pressure
45   USE sbcisf         ! surface boundary condition: ice shelf
46   USE sbcfwb         ! surface boundary condition: freshwater budget
47   USE icbstp         ! Icebergs
48   USE icb_oce  , ONLY : ln_passive_mode      ! iceberg interaction mode
49   USE traqsr         ! active tracers: light penetration
50   USE sbcwave        ! Wave module
51   USE bdy_oce   , ONLY: ln_bdy
52   USE usrdef_sbc     ! user defined: surface boundary condition
53   USE closea         ! closed sea
54   !
55   USE prtctl         ! Print control                    (prt_ctl routine)
56   USE iom            ! IOM library
57   USE in_out_manager ! I/O manager
58   USE lib_mpp        ! MPP library
59   USE timing         ! Timing
60   USE wet_dry
61   USE diurnal_bulk, ONLY:   ln_diurnal_only   ! diurnal SST diagnostic
62
63   IMPLICIT NONE
64   PRIVATE
65
66   PUBLIC   sbc        ! routine called by step.F90
67   PUBLIC   sbc_init   ! routine called by opa.F90
68
69   INTEGER ::   nsbc   ! type of surface boundary condition (deduced from namsbc informations)
70
71   !!----------------------------------------------------------------------
72   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
73   !! $Id$
74   !! Software governed by the CeCILL license (see ./LICENSE)
75   !!----------------------------------------------------------------------
76CONTAINS
77
78   SUBROUTINE sbc_init
79      !!---------------------------------------------------------------------
80      !!                    ***  ROUTINE sbc_init ***
81      !!
82      !! ** Purpose :   Initialisation of the ocean surface boundary computation
83      !!
84      !! ** Method  :   Read the namsbc namelist and set derived parameters
85      !!                Call init routines for all other SBC modules that have one
86      !!
87      !! ** Action  : - read namsbc parameters
88      !!              - nsbc: type of sbc
89      !!----------------------------------------------------------------------
90      INTEGER ::   ios, icpt                         ! local integer
91      LOGICAL ::   ll_purecpl, ll_opa, ll_not_nemo   ! local logical
92      !!
93      NAMELIST/namsbc/ nn_fsbc  ,                                                    &
94         &             ln_usr   , ln_flx   , ln_blk       ,                          &
95         &             ln_cpl   , ln_mixcpl, nn_components,                          &
96         &             nn_ice   , ln_ice_embd,                                       &
97         &             ln_traqsr, ln_dm2dc ,                                         &
98         &             ln_rnf   , nn_fwb   , ln_ssr   , ln_isf    , ln_apr_dyn ,     &
99         &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauwoc  , ln_stcor   ,     &
100         &             ln_tauw  , nn_lsm, nn_sdrift
101      !!----------------------------------------------------------------------
102      !
103      IF(lwp) THEN
104         WRITE(numout,*)
105         WRITE(numout,*) 'sbc_init : surface boundary condition setting'
106         WRITE(numout,*) '~~~~~~~~ '
107      ENDIF
108      !
109      !                       !**  read Surface Module namelist
110      REWIND( numnam_ref )          !* Namelist namsbc in reference namelist : Surface boundary
111      READ  ( numnam_ref, namsbc, IOSTAT = ios, ERR = 901)
112901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc in reference namelist', lwp )
113      REWIND( numnam_cfg )          !* Namelist namsbc in configuration namelist : Parameters of the run
114      READ  ( numnam_cfg, namsbc, IOSTAT = ios, ERR = 902 )
115902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc in configuration namelist', lwp )
116      IF(lwm) WRITE( numond, namsbc )
117      !
118#if defined key_mpp_mpi
119      ncom_fsbc = nn_fsbc    ! make nn_fsbc available for lib_mpp
120#endif
121      !                             !* overwrite namelist parameter using CPP key information
122#if defined key_agrif
123      IF( Agrif_Root() ) THEN                ! AGRIF zoom (cf r1242: possibility to run without ice in fine grid)
124         IF( lk_si3  )   nn_ice      = 2
125         IF( lk_cice )   nn_ice      = 3
126      ENDIF
127#else
128      IF( lk_si3  )   nn_ice      = 2
129      IF( lk_cice )   nn_ice      = 3
130#endif
131      !
132      IF(lwp) THEN                  !* Control print
133         WRITE(numout,*) '   Namelist namsbc (partly overwritten with CPP key setting)'
134         WRITE(numout,*) '      frequency update of sbc (and ice)             nn_fsbc       = ', nn_fsbc
135         WRITE(numout,*) '      Type of air-sea fluxes : '
136         WRITE(numout,*) '         user defined formulation                   ln_usr        = ', ln_usr
137         WRITE(numout,*) '         flux         formulation                   ln_flx        = ', ln_flx
138         WRITE(numout,*) '         bulk         formulation                   ln_blk        = ', ln_blk
139         WRITE(numout,*) '      Type of coupling (Ocean/Ice/Atmosphere) : '
140         WRITE(numout,*) '         ocean-atmosphere coupled formulation       ln_cpl        = ', ln_cpl
141         WRITE(numout,*) '         mixed forced-coupled     formulation       ln_mixcpl     = ', ln_mixcpl
142!!gm  lk_oasis is controlled by key_oasis3  ===>>>  It shoud be removed from the namelist
143         WRITE(numout,*) '         OASIS coupling (with atm or sas)           lk_oasis      = ', lk_oasis
144         WRITE(numout,*) '         components of your executable              nn_components = ', nn_components
145         WRITE(numout,*) '      Sea-ice : '
146         WRITE(numout,*) '         ice management in the sbc (=0/1/2/3)       nn_ice        = ', nn_ice
147         WRITE(numout,*) '         ice embedded into ocean                    ln_ice_embd   = ', ln_ice_embd
148         WRITE(numout,*) '      Misc. options of sbc : '
149         WRITE(numout,*) '         Light penetration in temperature Eq.       ln_traqsr     = ', ln_traqsr
150         WRITE(numout,*) '            daily mean to diurnal cycle qsr            ln_dm2dc   = ', ln_dm2dc
151         WRITE(numout,*) '         Sea Surface Restoring on SST and/or SSS    ln_ssr        = ', ln_ssr
152         WRITE(numout,*) '         FreshWater Budget control  (=0/1/2)        nn_fwb        = ', nn_fwb
153         WRITE(numout,*) '         Patm gradient added in ocean & ice Eqs.    ln_apr_dyn    = ', ln_apr_dyn
154         WRITE(numout,*) '         runoff / runoff mouths                     ln_rnf        = ', ln_rnf
155         WRITE(numout,*) '         iceshelf formulation                       ln_isf        = ', ln_isf
156         WRITE(numout,*) '         nb of iterations if land-sea-mask applied  nn_lsm        = ', nn_lsm
157         WRITE(numout,*) '         surface wave                               ln_wave       = ', ln_wave
158         WRITE(numout,*) '               Stokes drift corr. to vert. velocity ln_sdw        = ', ln_sdw
159         WRITE(numout,*) '                  vertical parametrization          nn_sdrift     = ', nn_sdrift
160         WRITE(numout,*) '               wave modified ocean stress           ln_tauwoc     = ', ln_tauwoc
161         WRITE(numout,*) '               wave modified ocean stress component ln_tauw       = ', ln_tauw
162         WRITE(numout,*) '               Stokes coriolis term                 ln_stcor      = ', ln_stcor
163         WRITE(numout,*) '               neutral drag coefficient (CORE,NCAR) ln_cdgw       = ', ln_cdgw
164      ENDIF
165      !
166      IF( .NOT.ln_wave ) THEN
167         ln_sdw = .false. ; ln_cdgw = .false. ; ln_tauwoc = .false. ; ln_tauw = .false. ; ln_stcor = .false.
168      ENDIF
169      IF( ln_sdw ) THEN
170         IF( .NOT.(nn_sdrift==jp_breivik_2014 .OR. nn_sdrift==jp_li_2017 .OR. nn_sdrift==jp_peakfr) ) &
171            CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' )
172      ENDIF
173      ll_st_bv2014  = ( nn_sdrift==jp_breivik_2014 )
174      ll_st_li2017  = ( nn_sdrift==jp_li_2017 )
175      ll_st_bv_li   = ( ll_st_bv2014 .OR. ll_st_li2017 )
176      ll_st_peakfr  = ( nn_sdrift==jp_peakfr )
177      IF( ln_tauwoc .AND. ln_tauw ) &
178         CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', &
179                                  '(ln_tauwoc=.true. and ln_tauw=.true.)' )
180      IF( ln_tauwoc ) &
181         CALL ctl_warn( 'You are subtracting the wave stress to the ocean (ln_tauwoc=.true.)' )
182      IF( ln_tauw ) &
183         CALL ctl_warn( 'The wave modified ocean stress components are used (ln_tauw=.true.) ', &
184                              'This will override any other specification of the ocean stress' )
185      !
186      IF( .NOT.ln_usr ) THEN     ! the model calendar needs some specificities (except in user defined case)
187         IF( MOD( rday , rdt ) /= 0. )   CALL ctl_stop( 'the time step must devide the number of second of in a day' )
188         IF( MOD( rday , 2.  ) /= 0. )   CALL ctl_stop( 'the number of second of in a day must be an even number'    )
189         IF( MOD( rdt  , 2.  ) /= 0. )   CALL ctl_stop( 'the time step (in second) must be an even number'           )
190      ENDIF
191      !                       !**  check option consistency
192      !
193      IF(lwp) WRITE(numout,*)       !* Single / Multi - executable (NEMO / OPA+SAS)
194      SELECT CASE( nn_components )
195      CASE( jp_iam_nemo )
196         IF(lwp) WRITE(numout,*) '   ==>>>   NEMO configured as a single executable (i.e. including both OPA and Surface module)'
197      CASE( jp_iam_opa  )
198         IF(lwp) WRITE(numout,*) '   ==>>>   Multi executable configuration. Here, OPA component'
199         IF( .NOT.lk_oasis )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' )
200         IF( ln_cpl        )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA'   )
201         IF( ln_mixcpl     )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' )
202      CASE( jp_iam_sas  )
203         IF(lwp) WRITE(numout,*) '   ==>>>   Multi executable configuration. Here, SAS component'
204         IF( .NOT.lk_oasis )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' )
205         IF( ln_mixcpl     )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' )
206      CASE DEFAULT
207         CALL ctl_stop( 'sbc_init : unsupported value for nn_components' )
208      END SELECT
209      !                             !* coupled options
210      IF( ln_cpl ) THEN
211         IF( .NOT. lk_oasis )   CALL ctl_stop( 'sbc_init : coupled mode with an atmosphere model (ln_cpl=T)',   &
212            &                                  '           required to defined key_oasis3' )
213      ENDIF
214      IF( ln_mixcpl ) THEN
215         IF( .NOT. lk_oasis )   CALL ctl_stop( 'sbc_init : mixed forced-coupled mode (ln_mixcpl=T) ',   &
216            &                                  '           required to defined key_oasis3' )
217         IF( .NOT.ln_cpl    )   CALL ctl_stop( 'sbc_init : mixed forced-coupled mode (ln_mixcpl=T) requires ln_cpl = T' )
218         IF( nn_components /= jp_iam_nemo )    &
219            &                   CALL ctl_stop( 'sbc_init : the mixed forced-coupled mode (ln_mixcpl=T) ',   &
220            &                                   '          not yet working with sas-opa coupling via oasis' )
221      ENDIF
222      !                             !* sea-ice
223      SELECT CASE( nn_ice )
224      CASE( 0 )                        !- no ice in the domain
225      CASE( 1 )                        !- Ice-cover climatology ("Ice-if" model) 
226      CASE( 2 )                        !- SI3  ice model
227      CASE( 3 )                        !- CICE ice model
228         IF( .NOT.( ln_blk .OR. ln_cpl ) )   CALL ctl_stop( 'sbc_init : CICE sea-ice model requires ln_blk or ln_cpl = T' )
229         IF( lk_agrif                    )   CALL ctl_stop( 'sbc_init : CICE sea-ice model not currently available with AGRIF' ) 
230      CASE DEFAULT                     !- not supported
231      END SELECT
232      !
233      !                       !**  allocate and set required variables
234      !
235      !                             !* allocate sbc arrays
236      IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'sbc_init : unable to allocate sbc_oce arrays' )
237#if ! defined key_si3 && ! defined key_cice
238      IF( sbc_ice_alloc() /= 0 )   CALL ctl_stop( 'sbc_init : unable to allocate sbc_ice arrays' )
239#endif
240      !
241      IF( .NOT.ln_isf ) THEN        !* No ice-shelf in the domain : allocate and set to zero
242         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' )
243         fwfisf  (:,:)   = 0._wp   ;   risf_tsc  (:,:,:) = 0._wp
244         fwfisf_b(:,:)   = 0._wp   ;   risf_tsc_b(:,:,:) = 0._wp
245      END IF
246      IF( nn_ice == 0 ) THEN        !* No sea-ice in the domain : ice fraction is always zero
247         IF( nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! except for OPA in SAS-OPA coupled case
248      ENDIF
249      !
250      sfx   (:,:) = 0._wp           !* salt flux due to freezing/melting
251      fmmflx(:,:) = 0._wp           !* freezing minus melting flux
252
253      taum(:,:) = 0._wp             !* wind stress module (needed in GLS in case of reduced restart)
254
255      !                          ! Choice of the Surface Boudary Condition (set nsbc)
256      IF( ln_dm2dc ) THEN           !* daily mean to diurnal cycle
257         nday_qsr = -1   ! allow initialization at the 1st call
258         IF( .NOT.( ln_flx .OR. ln_blk ) .AND. nn_components /= jp_iam_opa )   &
259            &   CALL ctl_stop( 'qsr diurnal cycle from daily values requires a flux or bulk formulation' )
260      ENDIF
261      !                             !* Choice of the Surface Boudary Condition
262      !                             (set nsbc)
263      !
264      ll_purecpl  = ln_cpl .AND. .NOT.ln_mixcpl
265      ll_opa      = nn_components == jp_iam_opa
266      ll_not_nemo = nn_components /= jp_iam_nemo
267      icpt = 0
268      !
269      IF( ln_usr          ) THEN   ;   nsbc = jp_usr     ; icpt = icpt + 1   ;   ENDIF       ! user defined         formulation
270      IF( ln_flx          ) THEN   ;   nsbc = jp_flx     ; icpt = icpt + 1   ;   ENDIF       ! flux                 formulation
271      IF( ln_blk          ) THEN   ;   nsbc = jp_blk     ; icpt = icpt + 1   ;   ENDIF       ! bulk                 formulation
272      IF( ll_purecpl      ) THEN   ;   nsbc = jp_purecpl ; icpt = icpt + 1   ;   ENDIF       ! Pure Coupled         formulation
273      IF( ll_opa          ) THEN   ;   nsbc = jp_none    ; icpt = icpt + 1   ;   ENDIF       ! opa coupling via SAS module
274      !
275      IF( icpt /= 1 )    CALL ctl_stop( 'sbc_init : choose ONE and only ONE sbc option' )
276      !
277      IF(lwp) THEN                     !- print the choice of surface flux formulation
278         WRITE(numout,*)
279         SELECT CASE( nsbc )
280         CASE( jp_usr     )   ;   WRITE(numout,*) '   ==>>>   user defined forcing formulation'
281         CASE( jp_flx     )   ;   WRITE(numout,*) '   ==>>>   flux formulation'
282         CASE( jp_blk     )   ;   WRITE(numout,*) '   ==>>>   bulk formulation'
283         CASE( jp_purecpl )   ;   WRITE(numout,*) '   ==>>>   pure coupled formulation'
284!!gm abusive use of jp_none ??   ===>>> need to be check and changed by adding a jp_sas parameter
285         CASE( jp_none    )   ;   WRITE(numout,*) '   ==>>>   OPA coupled to SAS via oasis'
286            IF( ln_mixcpl )       WRITE(numout,*) '               + forced-coupled mixed formulation'
287         END SELECT
288         IF( ll_not_nemo  )       WRITE(numout,*) '               + OASIS coupled SAS'
289      ENDIF
290      !
291      !                             !* OASIS initialization
292      !
293      IF( lk_oasis )   CALL sbc_cpl_init( nn_ice )   ! Must be done before: (1) first time step
294      !                                              !                      (2) the use of nn_fsbc
295      !     nn_fsbc initialization if OPA-SAS coupling via OASIS
296      !     SAS time-step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly
297      IF( nn_components /= jp_iam_nemo ) THEN
298         IF( nn_components == jp_iam_opa )   nn_fsbc = cpl_freq('O_SFLX') / NINT(rdt)
299         IF( nn_components == jp_iam_sas )   nn_fsbc = cpl_freq('I_SFLX') / NINT(rdt)
300         !
301         IF(lwp)THEN
302            WRITE(numout,*)
303            WRITE(numout,*)"   OPA-SAS coupled via OASIS : nn_fsbc re-defined from OASIS namcouple ", nn_fsbc
304            WRITE(numout,*)
305         ENDIF
306      ENDIF
307      !
308      !                             !* check consistency between model timeline and nn_fsbc
309      IF( MOD( nitend - nit000 + 1, nn_fsbc) /= 0 .OR.   &
310          MOD( nstock             , nn_fsbc) /= 0 ) THEN
311         WRITE(ctmp1,*) 'sbc_init : experiment length (', nitend - nit000 + 1, ') or nstock (', nstock,   &
312            &           ' is NOT a multiple of nn_fsbc (', nn_fsbc, ')'
313         CALL ctl_stop( ctmp1, 'Impossible to properly do model restart' )
314      ENDIF
315      !
316      IF( MOD( rday, REAL(nn_fsbc, wp) * rdt ) /= 0 )   &
317         &  CALL ctl_warn( 'sbc_init : nn_fsbc is NOT a multiple of the number of time steps in a day' )
318      !
319      IF( ln_dm2dc .AND. NINT(rday) / ( nn_fsbc * NINT(rdt) ) < 8  )   &
320         &   CALL ctl_warn( 'sbc_init : diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
321      !
322   
323      !                       !**  associated modules : initialization
324      !
325                          CALL sbc_ssm_init            ! Sea-surface mean fields initialization
326      !
327      IF( ln_blk      )   CALL sbc_blk_init            ! bulk formulae initialization
328
329      IF( ln_ssr      )   CALL sbc_ssr_init            ! Sea-Surface Restoring initialization
330      !
331      IF( ln_isf      )   CALL sbc_isf_init            ! Compute iceshelves
332      !
333                          CALL sbc_rnf_init            ! Runof initialization
334      !
335      IF( ln_apr_dyn )    CALL sbc_apr_init            ! Atmo Pressure Forcing initialization
336      !
337#if defined key_si3
338      IF( lk_agrif .AND. nn_ice == 0 ) THEN            ! allocate ice arrays in case agrif + ice-model + no-ice in child grid
339                          IF( sbc_ice_alloc() /= 0 )   CALL ctl_stop('STOP', 'sbc_ice_alloc : unable to allocate arrays' )
340      ELSEIF( nn_ice == 2 ) THEN
341                          CALL ice_init                ! ICE initialization
342      ENDIF
343#endif
344      IF( nn_ice == 3 )   CALL cice_sbc_init( nsbc )   ! CICE initialization
345      !
346      IF( ln_wave     )   CALL sbc_wave_init           ! surface wave initialisation
347      !
348      IF( lwxios ) THEN
349         CALL iom_set_rstw_var_active('utau_b')
350         CALL iom_set_rstw_var_active('vtau_b')
351         CALL iom_set_rstw_var_active('qns_b')
352         ! The 3D heat content due to qsr forcing is treated in traqsr
353         ! CALL iom_set_rstw_var_active('qsr_b')
354         CALL iom_set_rstw_var_active('emp_b')
355         CALL iom_set_rstw_var_active('sfx_b')
356      ENDIF
357
358   END SUBROUTINE sbc_init
359
360
361   SUBROUTINE sbc( kt )
362      !!---------------------------------------------------------------------
363      !!                    ***  ROUTINE sbc  ***
364      !!
365      !! ** Purpose :   provide at each time-step the ocean surface boundary
366      !!                condition (momentum, heat and freshwater fluxes)
367      !!
368      !! ** Method  :   blah blah  to be written ?????????
369      !!                CAUTION : never mask the surface stress field (tke sbc)
370      !!
371      !! ** Action  : - set the ocean surface boundary condition at before and now
372      !!                time step, i.e.
373      !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, sfx_b, qrp_b, erp_b
374      !!                utau  , vtau  , qns  , qsr  , emp  , sfx  , qrp  , erp
375      !!              - updte the ice fraction : fr_i
376      !!----------------------------------------------------------------------
377      INTEGER, INTENT(in) ::   kt   ! ocean time step
378      !
379      LOGICAL ::   ll_sas, ll_opa   ! local logical
380      !
381      REAL(wp) ::     zthscl        ! wd  tanh scale
382      REAL(wp), DIMENSION(jpi,jpj) ::  zwdht, zwght  ! wd dep over wd limit, wgt 
383
384      !!---------------------------------------------------------------------
385      !
386      IF( ln_timing )   CALL timing_start('sbc')
387      !
388      !                                            ! ---------------------------------------- !
389      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
390         !                                         ! ---------------------------------------- !
391         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields
392         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
393         qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine)
394         emp_b (:,:) = emp (:,:)
395         sfx_b (:,:) = sfx (:,:)
396         IF ( ln_rnf ) THEN
397            rnf_b    (:,:  ) = rnf    (:,:  )
398            rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
399         ENDIF
400         IF( ln_isf )  THEN
401            fwfisf_b  (:,:  ) = fwfisf  (:,:  )               
402            risf_tsc_b(:,:,:) = risf_tsc(:,:,:)             
403         ENDIF
404        !
405      ENDIF
406      !                                            ! ---------------------------------------- !
407      !                                            !        forcing field computation         !
408      !                                            ! ---------------------------------------- !
409      !
410      ll_sas = nn_components == jp_iam_sas               ! component flags
411      ll_opa = nn_components == jp_iam_opa
412      !
413      IF( .NOT.ll_sas )   CALL sbc_ssm ( kt )            ! mean ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
414      IF( ln_wave     )   CALL sbc_wave( kt )            ! surface waves
415
416      !
417      !                                            !==  sbc formulation  ==!
418      !                                                   
419      SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
420      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx)
421      CASE( jp_usr   )     ;   CALL usrdef_sbc_oce( kt )                    ! user defined formulation
422      CASE( jp_flx     )   ;   CALL sbc_flx       ( kt )                    ! flux formulation
423      CASE( jp_blk     )
424         IF( ll_sas    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: SAS receiving fields from OPA
425                               CALL sbc_blk       ( kt )                    ! bulk formulation for the ocean
426                               !
427      CASE( jp_purecpl )   ;   CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! pure coupled formulation
428      CASE( jp_none    )
429         IF( ll_opa    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: OPA receiving fields from SAS
430      END SELECT
431      !
432      IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing
433      !
434      IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )      ! Wind stress provided by waves
435      !
436      !                                            !==  Misc. Options  ==!
437      !
438      SELECT CASE( nn_ice )                                       ! Update heat and freshwater fluxes over sea-ice areas
439      CASE(  1 )   ;         CALL sbc_ice_if   ( kt )             ! Ice-cover climatology ("Ice-if" model)
440#if defined key_si3
441      CASE(  2 )   ;         CALL ice_stp  ( kt, nsbc )           ! SI3 ice model
442#endif
443      CASE(  3 )   ;         CALL sbc_ice_cice ( kt, nsbc )       ! CICE ice model
444      END SELECT
445
446      IF( ln_icebergs    )   THEN
447                                     CALL icb_stp( kt )           ! compute icebergs
448         ! icebergs may advect into haloes during the icb step and alter emp.
449         ! A lbc_lnk is necessary here to ensure restartability (#2113)
450         IF( .NOT. ln_passive_mode ) CALL lbc_lnk( 'sbcmod', emp, 'T', 1. ) ! ensure restartability with icebergs
451      ENDIF
452
453      IF( ln_isf         )   CALL sbc_isf( kt )                   ! compute iceshelves
454
455      IF( ln_rnf         )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes
456
457      IF( ln_ssr         )   CALL sbc_ssr( kt )                   ! add SST/SSS damping term
458
459      IF( nn_fwb    /= 0 )   CALL sbc_fwb( kt, nn_fwb, nn_fsbc )  ! control the freshwater budget
460
461      ! Special treatment of freshwater fluxes over closed seas in the model domain
462      ! Should not be run if ln_diurnal_only
463      IF( l_sbc_clo .AND. (.NOT. ln_diurnal_only) )   CALL sbc_clo( kt )   
464
465!!$!RBbug do not understand why see ticket 667
466!!$!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why.
467!!$      CALL lbc_lnk( 'sbcmod', emp, 'T', 1. )
468   IF ( ll_wd ) THEN     ! If near WAD point limit the flux for now
469    zthscl = atanh(rn_wd_sbcfra)                     ! taper frac default is .999
470    zwdht(:,:) = sshn(:,:) + ht_0(:,:) - rn_wdmin1   ! do this calc of water
471                                                     ! depth above wd limit once
472    WHERE( zwdht(:,:) <= 0.0 )
473            taum(:,:) = 0.0
474            utau(:,:) = 0.0
475            vtau(:,:) = 0.0
476            qns (:,:) = 0.0
477            qsr (:,:) = 0.0
478            emp (:,:) = min(emp(:,:),0.0) !can allow puddles to grow but not shrink
479            sfx (:,:) = 0.0
480    END WHERE
481    zwght(:,:) = tanh(zthscl*zwdht(:,:))
482    WHERE( zwdht(:,:) > 0.0  .and. zwdht(:,:) < rn_wd_sbcdep ) !  5 m hard limit here is arbitrary
483            qsr  (:,:) =  qsr(:,:)  * zwght(:,:)
484            qns  (:,:) =  qns(:,:)  * zwght(:,:)
485            taum (:,:) =  taum(:,:) * zwght(:,:)
486            utau (:,:) =  utau(:,:) * zwght(:,:)
487            vtau (:,:) =  vtau(:,:) * zwght(:,:)
488            sfx (:,:)  =  sfx(:,:)  * zwght(:,:)
489            emp  (:,:) =  emp(:,:)  * zwght(:,:)
490    END WHERE
491   ENDIF
492      !
493      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
494         !                                             ! ---------------------------------------- !
495         IF( ln_rstart .AND.    &                               !* Restart: read in restart file
496            & iom_varid( numror, 'utau_b', ldstop = .FALSE. ) > 0 ) THEN
497            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file'
498            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b, ldxios = lrxios )   ! before i-stress  (U-point)
499            CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b, ldxios = lrxios )   ! before j-stress  (V-point)
500            CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b, ldxios = lrxios  )   ! before non solar heat flux (T-point)
501            ! The 3D heat content due to qsr forcing is treated in traqsr
502            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b, ldxios = lrxios  ) ! before     solar heat flux (T-point)
503            CALL iom_get( numror, jpdom_autoglo, 'emp_b', emp_b, ldxios = lrxios  )    ! before     freshwater flux (T-point)
504            ! To ensure restart capability with 3.3x/3.4 restart files    !! to be removed in v3.6
505            IF( iom_varid( numror, 'sfx_b', ldstop = .FALSE. ) > 0 ) THEN
506               CALL iom_get( numror, jpdom_autoglo, 'sfx_b', sfx_b, ldxios = lrxios )  ! before salt flux (T-point)
507            ELSE
508               sfx_b (:,:) = sfx(:,:)
509            ENDIF
510         ELSE                                                   !* no restart: set from nit000 values
511            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
512            utau_b(:,:) = utau(:,:)
513            vtau_b(:,:) = vtau(:,:)
514            qns_b (:,:) = qns (:,:)
515            emp_b (:,:) = emp (:,:)
516            sfx_b (:,:) = sfx (:,:)
517         ENDIF
518      ENDIF
519      !                                                ! ---------------------------------------- !
520      IF( lrst_oce ) THEN                              !      Write in the ocean restart file     !
521         !                                             ! ---------------------------------------- !
522         IF(lwp) WRITE(numout,*)
523         IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   &
524            &                    'at it= ', kt,' date= ', ndastp
525         IF(lwp) WRITE(numout,*) '~~~~'
526         IF( lwxios ) CALL iom_swap(      cwxios_context          )
527         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau, ldxios = lwxios )
528         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau, ldxios = lwxios )
529         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns, ldxios = lwxios  )
530         ! The 3D heat content due to qsr forcing is treated in traqsr
531         ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
532         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp, ldxios = lwxios  )
533         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx, ldxios = lwxios  )
534         IF( lwxios ) CALL iom_swap(      cxios_context          )
535      ENDIF
536      !                                                ! ---------------------------------------- !
537      !                                                !        Outputs and control print         !
538      !                                                ! ---------------------------------------- !
539      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
540         CALL iom_put( "empmr"  , emp    - rnf )                ! upward water flux
541         CALL iom_put( "empbmr" , emp_b  - rnf )                ! before upward water flux ( needed to recalculate the time evolution of ssh in offline )
542         CALL iom_put( "saltflx", sfx  )                        ! downward salt flux (includes virtual salt flux beneath ice in linear free surface case)
543         CALL iom_put( "fmmflx", fmmflx  )                      ! Freezing-melting water flux
544         CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux
545         CALL iom_put( "qns"   , qns        )                   ! solar heat flux
546         CALL iom_put( "qsr"   ,       qsr  )                   ! solar heat flux
547         IF( nn_ice > 0 .OR. ll_opa )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction
548         CALL iom_put( "taum"  , taum       )                   ! wind stress module
549         CALL iom_put( "wspd"  , wndm       )                   ! wind speed  module over free ocean or leads in presence of sea-ice
550      ENDIF
551      !
552      CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at each time step in sea-ice)
553      CALL iom_put( "vtau", vtau )   ! j-wind stress
554      !
555      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
556         CALL prt_ctl(tab2d_1=fr_i              , clinfo1=' fr_i    - : ', mask1=tmask )
557         CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf), clinfo1=' emp-rnf - : ', mask1=tmask )
558         CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf), clinfo1=' sfx-rnf - : ', mask1=tmask )
559         CALL prt_ctl(tab2d_1=qns              , clinfo1=' qns      - : ', mask1=tmask )
560         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr      - : ', mask1=tmask )
561         CALL prt_ctl(tab3d_1=tmask            , clinfo1=' tmask    - : ', mask1=tmask, kdim=jpk )
562         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' sst      - : ', mask1=tmask, kdim=1   )
563         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sss      - : ', mask1=tmask, kdim=1   )
564         CALL prt_ctl(tab2d_1=utau             , clinfo1=' utau     - : ', mask1=umask,                      &
565            &         tab2d_2=vtau             , clinfo2=' vtau     - : ', mask2=vmask )
566      ENDIF
567
568      IF( kt == nitend )   CALL sbc_final         ! Close down surface module if necessary
569      !
570      IF( ln_timing )   CALL timing_stop('sbc')
571      !
572   END SUBROUTINE sbc
573
574
575   SUBROUTINE sbc_final
576      !!---------------------------------------------------------------------
577      !!                    ***  ROUTINE sbc_final  ***
578      !!
579      !! ** Purpose :   Finalize CICE (if used)
580      !!---------------------------------------------------------------------
581      !
582      IF( nn_ice == 3 )   CALL cice_sbc_final
583      !
584   END SUBROUTINE sbc_final
585
586   !!======================================================================
587END MODULE sbcmod
Note: See TracBrowser for help on using the repository browser.