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

Last change on this file since 12377 was 12377, checked in by acc, 8 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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