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 branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90 @ 7672

Last change on this file since 7672 was 7672, checked in by jcastill, 7 years ago

Further changes to remove the UM dependency in case of a coupled run, so that ocean can also run in coupled mode with a wave model only: until now, if running in forced/coupled mode (ln_mixcpl), the program expected forcing files and coupling fields from the atmosphere, and they were 'merged' together using a coupling map; if the coupling map was not provided, the coupling fields overwrote the input fields even if they were not actually coupled and did not have any valid information. If we were coupling to a wave model and not an atmosphere model, the forcing fields read from file were being overwritten but the atmosphere coupling fields, which did not contain any valid information.

Now, it is possible to couple in real mixed mode, where the fields can independently be read either from forcing files or from coupling (ocean coupled to an atmosphere model, a wave model, or both), and the coupling mapping is only needed if a field if both provided via coupling and a forcing file.

File size: 33.2 KB
Line 
1MODULE sbcmod
2   !!======================================================================
3   !!                       ***  MODULE  sbcmod  ***
4   !! Surface module :  provide to the ocean its surface boundary condition
5   !!======================================================================
6   !! History :  3.0  ! 2006-07  (G. Madec)  Original code
7   !!            3.1  ! 2008-08  (S. Masson, A. Caubel, E. Maisonnave, G. Madec) coupled interface
8   !!            3.3  ! 2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
9   !!            3.3  ! 2010-10  (S. Masson)  add diurnal cycle
10   !!            3.3  ! 2010-09  (D. Storkey) add ice boundary conditions (BDY)
11   !!             -   ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
12   !!             -   ! 2010-10  (J. Chanut, C. Bricaud, G. Madec)  add the surface pressure forcing
13   !!            3.4  ! 2011-11  (C. Harris) CICE added as an option
14   !!            3.5  ! 2012-11  (A. Coward, G. Madec) Rethink of heat, mass and salt surface fluxes
15   !!            3.6  ! 2014-11  (P. Mathiot, C. Harris) add ice shelves melting                   
16   !!----------------------------------------------------------------------
17
18   !!----------------------------------------------------------------------
19   !!   sbc_init       : read namsbc namelist
20   !!   sbc            : surface ocean momentum, heat and freshwater boundary conditions
21   !!----------------------------------------------------------------------
22   USE oce              ! ocean dynamics and tracers
23   USE dom_oce          ! ocean space and time domain
24   USE phycst           ! physical constants
25   USE sbc_oce          ! Surface boundary condition: ocean fields
26   USE trc_oce          ! shared ocean-passive tracers variables
27   USE sbc_ice          ! Surface boundary condition: ice fields
28   USE sbcdcy           ! surface boundary condition: diurnal cycle
29   USE sbcssm           ! surface boundary condition: sea-surface mean variables
30   USE sbcapr           ! surface boundary condition: atmospheric pressure
31   USE sbcana           ! surface boundary condition: analytical formulation
32   USE sbcflx           ! surface boundary condition: flux formulation
33   USE sbcblk_clio      ! surface boundary condition: bulk formulation : CLIO
34   USE sbcblk_core      ! surface boundary condition: bulk formulation : CORE
35   USE sbcblk_mfs       ! surface boundary condition: bulk formulation : MFS
36   USE sbcice_if        ! surface boundary condition: ice-if sea-ice model
37   USE sbcice_lim       ! surface boundary condition: LIM 3.0 sea-ice model
38   USE sbcice_lim_2     ! surface boundary condition: LIM 2.0 sea-ice model
39   USE sbcice_cice      ! surface boundary condition: CICE    sea-ice model
40   USE sbccpl           ! surface boundary condition: coupled florulation
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 sbcisf           ! surface boundary condition: ice shelf
45   USE sbcfwb           ! surface boundary condition: freshwater budget
46   USE closea           ! closed sea
47   USE icbstp           ! Icebergs!
48
49   USE prtctl           ! Print control                    (prt_ctl routine)
50   USE iom              ! IOM library
51   USE in_out_manager   ! I/O manager
52   USE lib_mpp          ! MPP library
53   USE timing           ! Timing
54   USE sbcwave          ! Wave module
55   USE bdy_par          ! Require lk_bdy
56
57   IMPLICIT NONE
58   PRIVATE
59
60   PUBLIC   sbc        ! routine called by step.F90
61   PUBLIC   sbc_init   ! routine called by opa.F90
62   
63   INTEGER ::   nsbc   ! type of surface boundary condition (deduced from namsbc informations)
64     
65   !! * Substitutions
66#  include "domzgr_substitute.h90"
67   !!----------------------------------------------------------------------
68   !! NEMO/OPA 4.0 , NEMO-consortium (2011)
69   !! $Id$
70   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
71   !!----------------------------------------------------------------------
72CONTAINS
73
74   SUBROUTINE sbc_init
75      !!---------------------------------------------------------------------
76      !!                    ***  ROUTINE sbc_init ***
77      !!
78      !! ** Purpose :   Initialisation of the ocean surface boundary computation
79      !!
80      !! ** Method  :   Read the namsbc namelist and set derived parameters
81      !!                Call init routines for all other SBC modules that have one
82      !!
83      !! ** Action  : - read namsbc parameters
84      !!              - nsbc: type of sbc
85      !!----------------------------------------------------------------------
86      INTEGER ::   icpt   ! local integer
87      !!
88      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl,   &
89         &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   &
90         &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   &
91         &             ln_tauoc  , ln_stcor  , nn_lsm, nn_limflx , nn_components, ln_cpl  ,   &
92         &             ln_wavcpl
93      INTEGER  ::   ios
94      INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
95      !!----------------------------------------------------------------------
96
97      IF(lwp) THEN
98         WRITE(numout,*)
99         WRITE(numout,*) 'sbc_init : surface boundary condition setting'
100         WRITE(numout,*) '~~~~~~~~ '
101      ENDIF
102
103      REWIND( numnam_ref )              ! Namelist namsbc in reference namelist : Surface boundary
104      READ  ( numnam_ref, namsbc, IOSTAT = ios, ERR = 901)
105901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc in reference namelist', lwp )
106
107      REWIND( numnam_cfg )              ! Namelist namsbc in configuration namelist : Parameters of the run
108      READ  ( numnam_cfg, namsbc, IOSTAT = ios, ERR = 902 )
109902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc in configuration namelist', lwp )
110      IF(lwm) WRITE ( numond, namsbc )
111
112      !                          ! overwrite namelist parameter using CPP key information
113      IF( Agrif_Root() ) THEN                ! AGRIF zoom
114        IF( lk_lim2 )   nn_ice      = 2
115        IF( lk_lim3 )   nn_ice      = 3
116        IF( lk_cice )   nn_ice      = 4
117      ENDIF
118      IF( cp_cfg == 'gyre' ) THEN            ! GYRE configuration
119          ln_ana      = .TRUE.   
120          nn_ice      =   0
121      ENDIF
122
123      IF(lwp) THEN               ! Control print
124         WRITE(numout,*) '        Namelist namsbc (partly overwritten with CPP key setting)'
125         WRITE(numout,*) '           frequency update of sbc (and ice)             nn_fsbc     = ', nn_fsbc
126         WRITE(numout,*) '           Type of sbc : '
127         WRITE(numout,*) '              analytical formulation                     ln_ana      = ', ln_ana
128         WRITE(numout,*) '              flux       formulation                     ln_flx      = ', ln_flx
129         WRITE(numout,*) '              CLIO bulk  formulation                     ln_blk_clio = ', ln_blk_clio
130         WRITE(numout,*) '              CORE bulk  formulation                     ln_blk_core = ', ln_blk_core
131         WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs  = ', ln_blk_mfs
132         WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl      = ', ln_cpl
133         WRITE(numout,*) '              forced-coupled atm mixed formulation       ln_mixcpl   = ', ln_mixcpl
134         WRITE(numout,*) '              forced-coupled wav mixed formulation       ln_wavcpl   = ', ln_wavcpl
135         WRITE(numout,*) '              wave physics                               ln_wave     = ', ln_wave 
136         WRITE(numout,*) '                 Stokes drift corr. to vert. velocity    ln_sdw      = ', ln_sdw 
137         WRITE(numout,*) '                 wave modified ocean stress              ln_tauoc    = ', ln_tauoc 
138         WRITE(numout,*) '                 Stokes coriolis term                    ln_stcor    = ', ln_stcor 
139         WRITE(numout,*) '                 neutral drag coefficient (CORE, MFS)    ln_cdgw     = ', ln_cdgw
140         WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis    = ', lk_oasis
141         WRITE(numout,*) '              components of your executable              nn_components = ', nn_components
142         WRITE(numout,*) '              Multicategory heat flux formulation (LIM3) nn_limflx   = ', nn_limflx
143         WRITE(numout,*) '           Misc. options of sbc : '
144         WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn
145         WRITE(numout,*) '              ice management in the sbc (=0/1/2/3)       nn_ice      = ', nn_ice 
146         WRITE(numout,*) '              ice-ocean embedded/levitating (=0/1/2)     nn_ice_embd = ', nn_ice_embd
147         WRITE(numout,*) '              daily mean to diurnal cycle qsr            ln_dm2dc    = ', ln_dm2dc 
148         WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf      = ', ln_rnf
149         WRITE(numout,*) '              iceshelf formulation                       nn_isf      = ', nn_isf
150         WRITE(numout,*) '              Sea Surface Restoring on SST and/or SSS    ln_ssr      = ', ln_ssr
151         WRITE(numout,*) '              FreshWater Budget control  (=0/1/2)        nn_fwb      = ', nn_fwb
152         WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea   = ', nn_closea
153         WRITE(numout,*) '              n. of iterations if land-sea-mask applied  nn_lsm      = ', nn_lsm
154      ENDIF
155
156      ! LIM3 Multi-category heat flux formulation
157      SELECT CASE ( nn_limflx)
158      CASE ( -1 )
159         IF(lwp) WRITE(numout,*) '              Use of per-category fluxes (nn_limflx = -1) '
160      CASE ( 0  )
161         IF(lwp) WRITE(numout,*) '              Average per-category fluxes (nn_limflx = 0) ' 
162      CASE ( 1  )
163         IF(lwp) WRITE(numout,*) '              Average then redistribute per-category fluxes (nn_limflx = 1) '
164      CASE ( 2  )
165         IF(lwp) WRITE(numout,*) '              Redistribute a single flux over categories (nn_limflx = 2) '
166      END SELECT
167      !
168      IF ( nn_components /= jp_iam_nemo .AND. .NOT. lk_oasis )   &
169         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' )
170      IF ( nn_components == jp_iam_opa .AND. ln_cpl )   &
171         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' )
172      IF ( nn_components == jp_iam_opa .AND. ( ln_mixcpl .OR. ln_wavcpl) )   &
173         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl or ln_wavcpl = T in OPA' )
174      IF ( ln_cpl .AND. .NOT. lk_oasis )    &
175         &      CALL ctl_stop( 'STOP', 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' )
176      IF( ( ln_mixcpl .OR. ln_wavcpl ) .AND. .NOT. lk_oasis )    &
177         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl or ln_wavcpl) requires the cpp key key_oasis3' )
178      IF( ( ln_mixcpl .OR. ln_wavcpl ) .AND. .NOT. ln_cpl )    &
179         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl or ln_wavcpl) requires ln_cpl = T' )
180      IF( ( ln_mixcpl .OR. ln_wavcpl ) .AND. nn_components /= jp_iam_nemo )    &
181         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl or ln_wavcpl) is not yet working with sas-opa coupling via oasis' )
182
183      !                              ! allocate sbc arrays
184      IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_oce arrays' )
185
186      !                          ! Checks:
187      IF( nn_isf .EQ. 0 ) THEN                      ! variable initialisation if no ice shelf
188         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' )
189         fwfisf  (:,:)   = 0.0_wp ; fwfisf_b  (:,:)   = 0.0_wp
190         risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp
191         rdivisf       = 0.0_wp
192      END IF
193      IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0.e0 ! no ice in the domain, ice fraction is always zero
194
195      sfx(:,:) = 0.0_wp                            ! the salt flux due to freezing/melting will be computed (i.e. will be non-zero)
196                                                   ! only if sea-ice is present
197 
198      fmmflx(:,:) = 0.0_wp                        ! freezing-melting array initialisation
199     
200      taum(:,:) = 0.0_wp                           ! Initialise taum for use in gls in case of reduced restart
201
202      !                                            ! restartability   
203      IF( ( nn_ice == 2 .OR. nn_ice ==3 ) .AND. .NOT.( ln_blk_clio .OR. ln_blk_core .OR. ln_cpl ) )   &
204         &   CALL ctl_stop( 'LIM sea-ice model requires a bulk formulation or coupled configuration' )
205      IF( nn_ice == 4 .AND. .NOT.( ln_blk_core .OR. ln_cpl ) )   &
206         &   CALL ctl_stop( 'CICE sea-ice model requires ln_blk_core or ln_cpl' )
207      IF( nn_ice == 4 .AND. lk_agrif )   &
208         &   CALL ctl_stop( 'CICE sea-ice model not currently available with AGRIF' )
209      IF( ( nn_ice == 3 .OR. nn_ice == 4 ) .AND. nn_ice_embd == 0 )   &
210         &   CALL ctl_stop( 'LIM3 and CICE sea-ice models require nn_ice_embd = 1 or 2' )
211      IF( ( nn_ice /= 3 ) .AND. ( nn_limflx >= 0 ) )   &
212         &   WRITE(numout,*) 'The nn_limflx>=0 option has no effect if sea ice model is not LIM3'
213      IF( ( nn_ice == 3 ) .AND. ( ln_cpl ) .AND. ( ( nn_limflx == -1 ) .OR. ( nn_limflx == 1 ) ) )   &
214         &   CALL ctl_stop( 'The chosen nn_limflx for LIM3 in coupled mode must be 0 or 2' )
215      IF( ( nn_ice == 3 ) .AND. ( .NOT. ln_cpl ) .AND. ( nn_limflx == 2 ) )   &
216         &   CALL ctl_stop( 'The chosen nn_limflx for LIM3 in forced mode cannot be 2' )
217
218      IF( ln_dm2dc )   nday_qsr = -1   ! initialisation flag
219
220      IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk_core ) .AND. nn_components /= jp_iam_opa )   &
221         &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' )
222     
223      IF ( ln_wave ) THEN
224         !Activated wave module but neither drag nor stokes drift activated
225         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) )   THEN 
226             CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F') 
227         !drag coefficient read from wave model definable only with mfs bulk formulae and core
228         ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN       
229             CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
230         ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN 
231             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
232         ENDIF
233      ELSE
234         IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) & 
235            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
236            &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
237            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 & 
238            &                  'or ocean stress modification due to waves (ln_tauoc=T) ',      & 
239            &                  'or Stokes-Coriolis term (ln_stcori=T)'  )
240      ENDIF 
241      !                          ! Choice of the Surface Boudary Condition (set nsbc)
242      ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl .AND. .NOT. ln_wavcpl
243      !
244      icpt = 0
245      IF( ln_ana          ) THEN   ;   nsbc = jp_ana     ; icpt = icpt + 1   ;   ENDIF       ! analytical           formulation
246      IF( ln_flx          ) THEN   ;   nsbc = jp_flx     ; icpt = icpt + 1   ;   ENDIF       ! flux                 formulation
247      IF( ln_blk_clio     ) THEN   ;   nsbc = jp_clio    ; icpt = icpt + 1   ;   ENDIF       ! CLIO bulk            formulation
248      IF( ln_blk_core     ) THEN   ;   nsbc = jp_core    ; icpt = icpt + 1   ;   ENDIF       ! CORE bulk            formulation
249      IF( ln_blk_mfs      ) THEN   ;   nsbc = jp_mfs     ; icpt = icpt + 1   ;   ENDIF       ! MFS  bulk            formulation
250      IF( ll_purecpl      ) THEN   ;   nsbc = jp_purecpl ; icpt = icpt + 1   ;   ENDIF       ! Pure Coupled         formulation
251      IF( cp_cfg == 'gyre') THEN   ;   nsbc = jp_gyre                        ;   ENDIF       ! GYRE analytical      formulation
252      IF( nn_components == jp_iam_opa )   &
253         &                  THEN   ;   nsbc = jp_none    ; icpt = icpt + 1   ;   ENDIF       ! opa coupling via SAS module
254      IF( lk_esopa        )            nsbc = jp_esopa                                       ! esopa test, ALL formulations
255      !
256      IF( icpt /= 1 .AND. .NOT.lk_esopa ) THEN
257         WRITE(numout,*)
258         WRITE(numout,*) '           E R R O R in setting the sbc, one and only one namelist/CPP key option '
259         WRITE(numout,*) '                     must be choosen. You choose ', icpt, ' option(s)'
260         WRITE(numout,*) '                     We stop'
261         nstop = nstop + 1
262      ENDIF
263      IF(lwp) THEN
264         WRITE(numout,*)
265         IF( nsbc == jp_esopa   )   WRITE(numout,*) '              ESOPA test All surface boundary conditions'
266         IF( nsbc == jp_gyre    )   WRITE(numout,*) '              GYRE analytical formulation'
267         IF( nsbc == jp_ana     )   WRITE(numout,*) '              analytical formulation'
268         IF( nsbc == jp_flx     )   WRITE(numout,*) '              flux formulation'
269         IF( nsbc == jp_clio    )   WRITE(numout,*) '              CLIO bulk formulation'
270         IF( nsbc == jp_core    )   WRITE(numout,*) '              CORE bulk formulation'
271         IF( nsbc == jp_purecpl )   WRITE(numout,*) '              pure coupled formulation'
272         IF( nsbc == jp_mfs     )   WRITE(numout,*) '              MFS Bulk formulation'
273         IF( nsbc == jp_none    )   WRITE(numout,*) '              OPA coupled to SAS via oasis'
274         IF( ln_mixcpl          )   WRITE(numout,*) '              + forced-coupled mixed atm formulation'
275         IF( ln_wavcpl          )   WRITE(numout,*) '              + forced-coupled mixed wav formulation'
276         IF( nn_components/= jp_iam_nemo )  &
277            &                       WRITE(numout,*) '              + OASIS coupled SAS'
278      ENDIF
279      !
280      IF( lk_oasis )   CALL sbc_cpl_init (nn_ice)   ! OASIS initialisation. must be done before: (1) first time step
281      !                                                     !                                            (2) the use of nn_fsbc
282
283!     nn_fsbc initialization if OPA-SAS coupling via OASIS
284!     sas model time step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly
285      IF ( nn_components /= jp_iam_nemo ) THEN
286
287         IF ( nn_components == jp_iam_opa ) nn_fsbc = cpl_freq('O_SFLX') / NINT(rdt)
288         IF ( nn_components == jp_iam_sas ) nn_fsbc = cpl_freq('I_SFLX') / NINT(rdt)
289         !
290         IF(lwp)THEN
291            WRITE(numout,*)
292            WRITE(numout,*)"   OPA-SAS coupled via OASIS : nn_fsbc re-defined from OASIS namcouple ", nn_fsbc
293            WRITE(numout,*)
294         ENDIF
295      ENDIF
296
297      IF( MOD( nitend - nit000 + 1, nn_fsbc) /= 0 .OR.   &
298          MOD( nstock             , nn_fsbc) /= 0 ) THEN
299         WRITE(ctmp1,*) 'experiment length (', nitend - nit000 + 1, ') or nstock (', nstock,   &
300            &           ' is NOT a multiple of nn_fsbc (', nn_fsbc, ')'
301         CALL ctl_stop( ctmp1, 'Impossible to properly do model restart' )
302      ENDIF
303      !
304      IF( MOD( rday, REAL(nn_fsbc, wp) * rdt ) /= 0 )   &
305         &  CALL ctl_warn( 'nn_fsbc is NOT a multiple of the number of time steps in a day' )
306      !
307      IF( ln_dm2dc .AND. ( ( NINT(rday) / ( nn_fsbc * NINT(rdt) ) )  < 8 ) )   &
308         &   CALL ctl_warn( 'diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
309
310                               CALL sbc_ssm_init               ! Sea-surface mean fields initialisation
311      !
312      IF( ln_ssr           )   CALL sbc_ssr_init               ! Sea-Surface Restoring initialisation
313      !
314                               CALL sbc_rnf_init               ! Runof initialisation
315      !
316      IF( nn_ice == 3      )   CALL sbc_lim_init               ! LIM3 initialisation
317
318      IF( nn_ice == 4      )   CALL cice_sbc_init( nsbc )      ! CICE initialisation
319      !
320      IF( ln_wave          )   CALL sbc_wave_init              ! surface wave initialisation
321      !
322     
323   END SUBROUTINE sbc_init
324
325
326   SUBROUTINE sbc( kt )
327      !!---------------------------------------------------------------------
328      !!                    ***  ROUTINE sbc  ***
329      !!             
330      !! ** Purpose :   provide at each time-step the ocean surface boundary
331      !!                condition (momentum, heat and freshwater fluxes)
332      !!
333      !! ** Method  :   blah blah  to be written ?????????
334      !!                CAUTION : never mask the surface stress field (tke sbc)
335      !!
336      !! ** Action  : - set the ocean surface boundary condition at before and now
337      !!                time step, i.e. 
338      !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, sfx_b, qrp_b, erp_b
339      !!                utau  , vtau  , qns  , qsr  , emp  , sfx  , qrp  , erp
340      !!              - updte the ice fraction : fr_i
341      !!----------------------------------------------------------------------
342      INTEGER, INTENT(in) ::   kt       ! ocean time step
343      !!---------------------------------------------------------------------
344      !
345      IF( nn_timing == 1 )  CALL timing_start('sbc')
346      !
347      !                                            ! ---------------------------------------- !
348      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
349         !                                         ! ---------------------------------------- !
350         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields
351         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
352         qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine)
353         ! The 3D heat content due to qsr forcing is treated in traqsr
354         ! qsr_b (:,:) = qsr (:,:)
355         emp_b(:,:) = emp(:,:)
356         sfx_b(:,:) = sfx(:,:)
357      ENDIF
358      !                                            ! ---------------------------------------- !
359      !                                            !        forcing field computation         !
360      !                                            ! ---------------------------------------- !
361      !
362      IF ( .NOT. lk_bdy ) then
363         IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc
364      ENDIF
365                                                         ! (caution called before sbc_ssm)
366      !
367      IF( nn_components /= jp_iam_sas )   CALL sbc_ssm( kt )   ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
368      !                                                        ! averaged over nf_sbc time-step
369
370      IF (ln_wave) CALL sbc_wave( kt )
371                                                   !==  sbc formulation  ==!
372                                                           
373      SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
374      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx)
375      CASE( jp_gyre  )   ;   CALL sbc_gyre    ( kt )                    ! analytical formulation : GYRE configuration
376      CASE( jp_ana   )   ;   CALL sbc_ana     ( kt )                    ! analytical formulation : uniform sbc
377      CASE( jp_flx   )   ;   CALL sbc_flx     ( kt )                    ! flux formulation
378      CASE( jp_clio  )   ;   CALL sbc_blk_clio( kt )                    ! bulk formulation : CLIO for the ocean
379      CASE( jp_core  )   
380         IF( nn_components == jp_iam_sas ) &
381            &                CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: SAS receiving fields from OPA
382                             CALL sbc_blk_core( kt )                    ! bulk formulation : CORE for the ocean
383                                                                        ! from oce: sea surface variables (sst_m, sss_m,  ssu_m,  ssv_m)
384      CASE( jp_purecpl )  ;  CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! pure coupled formulation
385                                                                        !
386      CASE( jp_mfs   )   ;   CALL sbc_blk_mfs ( kt )                    ! bulk formulation : MFS for the ocean
387      CASE( jp_none  ) 
388         IF( nn_components == jp_iam_opa ) &
389                             CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: OPA receiving fields from SAS
390      CASE( jp_esopa )                               
391                             CALL sbc_ana     ( kt )                    ! ESOPA, test ALL the formulations
392                             CALL sbc_gyre    ( kt )                    !
393                             CALL sbc_flx     ( kt )                    !
394                             CALL sbc_blk_clio( kt )                    !
395                             CALL sbc_blk_core( kt )                    !
396                             CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   !
397      END SELECT
398
399      IF( ln_mixcpl .OR. ln_wavcpl )  CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing
400
401      IF ( ln_wave .AND. ln_tauoc) THEN                 ! Wave stress subctracted
402            utau(:,:) = utau(:,:)*tauoc_wave(:,:) 
403            vtau(:,:) = vtau(:,:)*tauoc_wave(:,:) 
404            taum(:,:) = taum(:,:)*tauoc_wave(:,:) 
405      !
406            SELECT CASE( nsbc ) 
407            CASE(  0,1,2,3,5,-1 )  ; 
408                IF(lwp .AND. kt == nit000 ) WRITE(numout,*) 'WARNING: You are subtracting the wave stress to the ocean. &
409                        & If not requested select ln_tauoc=.false' 
410            END SELECT 
411      !
412      END IF
413
414      !                                            !==  Misc. Options  ==!
415     
416      SELECT CASE( nn_ice )                                       ! Update heat and freshwater fluxes over sea-ice areas
417      CASE(  1 )   ;         CALL sbc_ice_if   ( kt )                ! Ice-cover climatology ("Ice-if" model)
418      CASE(  2 )   ;         CALL sbc_ice_lim_2( kt, nsbc )          ! LIM-2 ice model
419      CASE(  3 )   ;         CALL sbc_ice_lim  ( kt, nsbc )          ! LIM-3 ice model
420      CASE(  4 )   ;         CALL sbc_ice_cice ( kt, nsbc )          ! CICE ice model
421      END SELECT                                             
422
423      IF( ln_icebergs    )   CALL icb_stp( kt )                   ! compute icebergs
424
425      IF( nn_isf   /= 0  )   CALL sbc_isf( kt )                    ! compute iceshelves
426
427      IF( ln_rnf         )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes
428 
429      IF( ln_ssr         )   CALL sbc_ssr( kt )                   ! add SST/SSS damping term
430
431      IF( nn_fwb    /= 0 )   CALL sbc_fwb( kt, nn_fwb, nn_fsbc )  ! control the freshwater budget
432
433      IF( nn_closea == 1 )   CALL sbc_clo( kt )                   ! treatment of closed sea in the model domain
434      !                                                           ! (update freshwater fluxes)
435!RBbug do not understand why see ticket 667
436!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why.
437      CALL lbc_lnk( emp, 'T', 1. )
438      !
439      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
440         !                                             ! ---------------------------------------- !
441         IF( ln_rstart .AND.    &                               !* Restart: read in restart file
442            & iom_varid( numror, 'utau_b', ldstop = .FALSE. ) > 0 ) THEN
443            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file'
444            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b )   ! before i-stress  (U-point)
445            CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b )   ! before j-stress  (V-point)
446            CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b  )   ! before non solar heat flux (T-point)
447            ! The 3D heat content due to qsr forcing is treated in traqsr
448            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  ) ! before     solar heat flux (T-point)
449            CALL iom_get( numror, jpdom_autoglo, 'emp_b', emp_b  )    ! before     freshwater flux (T-point)
450            ! To ensure restart capability with 3.3x/3.4 restart files    !! to be removed in v3.6
451            IF( iom_varid( numror, 'sfx_b', ldstop = .FALSE. ) > 0 ) THEN
452               CALL iom_get( numror, jpdom_autoglo, 'sfx_b', sfx_b )  ! before salt flux (T-point)
453            ELSE
454               sfx_b (:,:) = sfx(:,:)
455            ENDIF
456         ELSE                                                   !* no restart: set from nit000 values
457            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
458            utau_b(:,:) = utau(:,:) 
459            vtau_b(:,:) = vtau(:,:)
460            qns_b (:,:) = qns (:,:)
461            emp_b (:,:) = emp(:,:)
462            sfx_b (:,:) = sfx(:,:)
463         ENDIF
464      ENDIF
465      !                                                ! ---------------------------------------- !
466      IF( lrst_oce ) THEN                              !      Write in the ocean restart file     !
467         !                                             ! ---------------------------------------- !
468         IF(lwp) WRITE(numout,*)
469         IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   &
470            &                    'at it= ', kt,' date= ', ndastp
471         IF(lwp) WRITE(numout,*) '~~~~'
472         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau )
473         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau )
474         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  )
475         ! The 3D heat content due to qsr forcing is treated in traqsr
476         ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
477         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  )
478         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx  )
479      ENDIF
480
481      !                                                ! ---------------------------------------- !
482      !                                                !        Outputs and control print         !
483      !                                                ! ---------------------------------------- !
484      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
485         CALL iom_put( "empmr" , emp  - rnf )                   ! upward water flux
486         CALL iom_put( "saltflx", sfx  )                        ! downward salt flux 
487                                                                ! (includes virtual salt flux beneath ice
488                                                                ! in linear free surface case)
489         CALL iom_put( "fmmflx", fmmflx  )                      ! Freezing-melting water flux
490         CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux
491         CALL iom_put( "qns"   , qns        )                   ! solar heat flux
492         CALL iom_put( "qsr"   ,       qsr  )                   ! solar heat flux
493         IF( nn_ice > 0 .OR. nn_components == jp_iam_opa )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction
494         CALL iom_put( "taum"  , taum       )                   ! wind stress module
495         CALL iom_put( "wspd"  , wndm       )                   ! wind speed  module over free ocean or leads in presence of sea-ice
496      ENDIF
497      !
498      CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at
499      CALL iom_put( "vtau", vtau )   ! j-wind stress    each time step in sea-ice)
500      !
501      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
502         CALL prt_ctl(tab2d_1=fr_i              , clinfo1=' fr_i     - : ', mask1=tmask, ovlap=1 )
503         CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf), clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 )
504         CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf), clinfo1=' sfx-rnf  - : ', mask1=tmask, ovlap=1 )
505         CALL prt_ctl(tab2d_1=qns              , clinfo1=' qns      - : ', mask1=tmask, ovlap=1 )
506         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr      - : ', mask1=tmask, ovlap=1 )
507         CALL prt_ctl(tab3d_1=tmask            , clinfo1=' tmask    - : ', mask1=tmask, ovlap=1, kdim=jpk )
508         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' sst      - : ', mask1=tmask, ovlap=1, kdim=1   )
509         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sss      - : ', mask1=tmask, ovlap=1, kdim=1   )
510         CALL prt_ctl(tab2d_1=utau             , clinfo1=' utau     - : ', mask1=umask,                      &
511            &         tab2d_2=vtau             , clinfo2=' vtau     - : ', mask2=vmask, ovlap=1 )
512      ENDIF
513
514      IF( kt == nitend )   CALL sbc_final         ! Close down surface module if necessary
515      !
516      IF( nn_timing == 1 )  CALL timing_stop('sbc')
517      !
518   END SUBROUTINE sbc
519
520
521   SUBROUTINE sbc_final
522      !!---------------------------------------------------------------------
523      !!                    ***  ROUTINE sbc_final  ***
524      !!
525      !! ** Purpose :   Finalize CICE (if used)
526      !!---------------------------------------------------------------------
527      !
528      IF( nn_ice == 4 )   CALL cice_sbc_final
529      !
530   END SUBROUTINE sbc_final
531
532   !!======================================================================
533END MODULE sbcmod
Note: See TracBrowser for help on using the repository browser.