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.
sbcmod.F90 in NEMO/trunk/src/OCE/SBC – NEMO

source: NEMO/trunk/src/OCE/SBC/sbcmod.F90 @ 9940

Last change on this file since 9940 was 9940, checked in by acc, 6 years ago

Changes to restore restartability and reproducibility with ORCA2_ICE_PISCES SETTE tests with active icebergs. Changes to icbthm.F90 are purely cosmetic or improvements to code efficiency. The key change is the reintroduction of an lbc_lnk for emp in sbcmod after icb_stp. Icebergs that advect into haloes during icb_stp can melt and alter emp (if ln_passive_mode is .false.). These changes are lost across restarts without the halo exchange. This solution is not ideal but no alterantive is apparent with the current icb algorithm. This change restores both restartability and reproducibility with a 990 timestep (495 timestep restart) set of tests. See #2113 for details

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