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/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcmod.F90 @ 10989

Last change on this file since 10989 was 10989, checked in by acc, 4 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert NST routines in preparation for getting AGRIF back up and running. AGRIF conv stage now works but requires some renaming of recently changes DIU modules (included in this commit). AGRIF compile and link stage not yet working (agrif routines need to be passed the time-level indices) but non-AGRIF SETTE tests are all OK

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