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/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90 @ 2797

Last change on this file since 2797 was 2797, checked in by davestorkey, 13 years ago

Delete BDY module and first implementation of new OBC module.

  1. Initial restructuring.
  2. Use fldread to read open boundary data.
  • Property svn:keywords set to Id
File size: 20.7 KB
Line 
1MODULE sbcmod
2   !!======================================================================
3   !!                       ***  MODULE  sbcmod  ***
4   !! Surface module :  provide to the ocean its surface boundary condition
5   !!======================================================================
6   !! History :  3.0  ! 2006-07  (G. Madec)  Original code
7   !!            3.1  ! 2008-08  (S. Masson, A. Caubel, E. Maisonnave, G. Madec) coupled interface
8   !!            3.3  ! 2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
9   !!            3.3  ! 2010-10  (S. Masson)  add diurnal cycle
10   !!            3.3  ! 2010-09  (D. Storkey) add ice boundary conditions (BDY)
11   !!             -   ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
12   !!             -   ! 2010-10  (J. Chanut, C. Bricaud, G. Madec)  add the surface pressure forcing
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   sbc_init       : read namsbc namelist
17   !!   sbc            : surface ocean momentum, heat and freshwater boundary conditions
18   !!----------------------------------------------------------------------
19   USE oce              ! ocean dynamics and tracers
20   USE dom_oce          ! ocean space and time domain
21   USE phycst           ! physical constants
22   USE sbc_oce          ! Surface boundary condition: ocean fields
23   USE sbc_ice          ! Surface boundary condition: ice fields
24   USE sbcdcy           ! surface boundary condition: diurnal cycle
25   USE sbcssm           ! surface boundary condition: sea-surface mean variables
26   USE sbcapr           ! surface boundary condition: atmospheric pressure
27   USE sbcana           ! surface boundary condition: analytical formulation
28   USE sbcflx           ! surface boundary condition: flux formulation
29   USE sbcblk_clio      ! surface boundary condition: bulk formulation : CLIO
30   USE sbcblk_core      ! surface boundary condition: bulk formulation : CORE
31   USE sbcice_if        ! surface boundary condition: ice-if sea-ice model
32   USE sbcice_lim       ! surface boundary condition: LIM 3.0 sea-ice model
33   USE sbcice_lim_2     ! surface boundary condition: LIM 2.0 sea-ice model
34   USE sbccpl           ! surface boundary condition: coupled florulation
35   USE cpl_oasis3, ONLY:lk_cpl      ! are we in coupled mode?
36   USE sbcssr           ! surface boundary condition: sea surface restoring
37   USE sbcrnf           ! surface boundary condition: runoffs
38   USE sbcfwb           ! surface boundary condition: freshwater budget
39   USE closea           ! closed sea
40   USE obc_par          ! for lk_obc
41   USE obcice_lim2      ! unstructured open boundary data  (obc_ice_lim_2 routine)
42
43   USE prtctl           ! Print control                    (prt_ctl routine)
44   USE restart          ! ocean restart
45   USE iom              ! IOM library
46   USE in_out_manager   ! I/O manager
47   USE lib_mpp          ! MPP library
48
49   IMPLICIT NONE
50   PRIVATE
51
52   PUBLIC   sbc        ! routine called by step.F90
53   PUBLIC   sbc_init   ! routine called by opa.F90
54   
55   INTEGER ::   nsbc   ! type of surface boundary condition (deduced from namsbc informations)
56     
57   !! * Substitutions
58#  include "domzgr_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/OPA 4.0 , NEMO-consortium (2011)
61   !! $Id$
62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE sbc_init
67      !!---------------------------------------------------------------------
68      !!                    ***  ROUTINE sbc_init ***
69      !!
70      !! ** Purpose :   Initialisation of the ocean surface boundary computation
71      !!
72      !! ** Method  :   Read the namsbc namelist and set derived parameters
73      !!
74      !! ** Action  : - read namsbc parameters
75      !!              - nsbc: type of sbc
76      !!----------------------------------------------------------------------
77      INTEGER ::   icpt   ! local integer
78      !!
79      NAMELIST/namsbc/ nn_fsbc   , ln_ana , ln_flx  , ln_blk_clio, ln_blk_core, ln_cpl,   &
80         &             ln_apr_dyn, nn_ice , ln_dm2dc, ln_rnf     , ln_ssr     , nn_fwb
81      !!----------------------------------------------------------------------
82
83      IF(lwp) THEN
84         WRITE(numout,*)
85         WRITE(numout,*) 'sbc_init : surface boundary condition setting'
86         WRITE(numout,*) '~~~~~~~~ '
87      ENDIF
88
89      REWIND( numnam )           ! Read Namelist namsbc
90      READ  ( numnam, namsbc )
91
92      !                          ! overwrite namelist parameter using CPP key information
93      IF( Agrif_Root() ) THEN                ! AGRIF zoom
94        IF( lk_lim2 )   nn_ice      = 2
95        IF( lk_lim3 )   nn_ice      = 3
96      ENDIF
97      IF( cp_cfg == 'gyre' ) THEN            ! GYRE configuration
98          ln_ana      = .TRUE.   
99          nn_ice      =   0
100      ENDIF
101     
102      IF(lwp) THEN               ! Control print
103         WRITE(numout,*) '        Namelist namsbc (partly overwritten with CPP key setting)'
104         WRITE(numout,*) '           frequency update of sbc (and ice)             nn_fsbc     = ', nn_fsbc
105         WRITE(numout,*) '           Type of sbc : '
106         WRITE(numout,*) '              analytical formulation                     ln_ana      = ', ln_ana
107         WRITE(numout,*) '              flux       formulation                     ln_flx      = ', ln_flx
108         WRITE(numout,*) '              CLIO bulk  formulation                     ln_blk_clio = ', ln_blk_clio
109         WRITE(numout,*) '              CLIO bulk  formulation                     ln_blk_core = ', ln_blk_core
110         WRITE(numout,*) '              coupled    formulation (T if key_sbc_cpl)  ln_cpl      = ', ln_cpl
111         WRITE(numout,*) '           Misc. options of sbc : '
112         WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn
113         WRITE(numout,*) '              ice management in the sbc (=0/1/2/3)       nn_ice      = ', nn_ice 
114         WRITE(numout,*) '              daily mean to diurnal cycle qsr            ln_dm2dc    = ', ln_dm2dc 
115         WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf      = ', ln_rnf
116         WRITE(numout,*) '              Sea Surface Restoring on SST and/or SSS    ln_ssr      = ', ln_ssr
117         WRITE(numout,*) '              FreshWater Budget control  (=0/1/2)        nn_fwb      = ', nn_fwb
118         WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea   = ', nn_closea
119      ENDIF
120
121      !                              ! allocate sbc arrays
122      IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_oce arrays' )
123
124      !                          ! Checks:
125      IF( .NOT. ln_rnf ) THEN                      ! no specific treatment in vicinity of river mouths
126         ln_rnf_mouth  = .false.                     
127         IF( sbc_rnf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_rnf arrays' )
128         nkrnf         = 0
129         rnf     (:,:) = 0.e0
130         rnfmsk  (:,:) = 0.e0
131         rnfmsk_z(:)   = 0.e0
132      ENDIF
133      IF( nn_ice == 0  )   fr_i(:,:) = 0.e0        ! no ice in the domain, ice fraction is always zero
134
135      !                                            ! restartability   
136      IF( MOD( nitend - nit000 + 1, nn_fsbc) /= 0 .OR.   &
137          MOD( nstock             , nn_fsbc) /= 0 ) THEN
138         WRITE(ctmp1,*) 'experiment length (', nitend - nit000 + 1, ') or nstock (', nstock,   &
139            &           ' is NOT a multiple of nn_fsbc (', nn_fsbc, ')'
140         CALL ctl_stop( ctmp1, 'Impossible to properly do model restart' )
141      ENDIF
142      !
143      IF( MOD( rday, REAL(nn_fsbc, wp) * rdt ) /= 0 )   &
144         &  CALL ctl_warn( 'nn_fsbc is NOT a multiple of the number of time steps in a day' )
145      !
146      IF( nn_ice == 2 .AND. .NOT.( ln_blk_clio .OR. ln_blk_core .OR. lk_cpl ) )   &
147         &   CALL ctl_stop( 'sea-ice model requires a bulk formulation or coupled configuration' )
148     
149      IF( ln_dm2dc )   nday_qsr = -1   ! initialisation flag
150
151      IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk_core ) )   &
152         &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' )
153     
154      IF( ln_dm2dc .AND. ( ( NINT(rday) / ( nn_fsbc * NINT(rdt) ) )  < 8 ) )   &
155         &   CALL ctl_warn( 'diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
156     
157      !                          ! Choice of the Surface Boudary Condition (set nsbc)
158      icpt = 0
159      IF( ln_ana          ) THEN   ;   nsbc =  1   ; icpt = icpt + 1   ;   ENDIF       ! analytical      formulation
160      IF( ln_flx          ) THEN   ;   nsbc =  2   ; icpt = icpt + 1   ;   ENDIF       ! flux            formulation
161      IF( ln_blk_clio     ) THEN   ;   nsbc =  3   ; icpt = icpt + 1   ;   ENDIF       ! CLIO bulk       formulation
162      IF( ln_blk_core     ) THEN   ;   nsbc =  4   ; icpt = icpt + 1   ;   ENDIF       ! CORE bulk       formulation
163      IF( ln_cpl          ) THEN   ;   nsbc =  5   ; icpt = icpt + 1   ;   ENDIF       ! Coupled         formulation
164      IF( cp_cfg == 'gyre') THEN   ;   nsbc =  0                       ;   ENDIF       ! GYRE analytical formulation
165      IF( lk_esopa        )            nsbc = -1                                       ! esopa test, ALL formulations
166      !
167      IF( icpt /= 1 .AND. .NOT.lk_esopa ) THEN
168         WRITE(numout,*)
169         WRITE(numout,*) '           E R R O R in setting the sbc, one and only one namelist/CPP key option '
170         WRITE(numout,*) '                     must be choosen. You choose ', icpt, ' option(s)'
171         WRITE(numout,*) '                     We stop'
172         nstop = nstop + 1
173      ENDIF
174      IF(lwp) THEN
175         WRITE(numout,*)
176         IF( nsbc == -1 )   WRITE(numout,*) '              ESOPA test All surface boundary conditions'
177         IF( nsbc ==  0 )   WRITE(numout,*) '              GYRE analytical formulation'
178         IF( nsbc ==  1 )   WRITE(numout,*) '              analytical formulation'
179         IF( nsbc ==  2 )   WRITE(numout,*) '              flux formulation'
180         IF( nsbc ==  3 )   WRITE(numout,*) '              CLIO bulk formulation'
181         IF( nsbc ==  4 )   WRITE(numout,*) '              CORE bulk formulation'
182         IF( nsbc ==  5 )   WRITE(numout,*) '              coupled formulation'
183      ENDIF
184      !
185   END SUBROUTINE sbc_init
186
187
188   SUBROUTINE sbc( kt )
189      !!---------------------------------------------------------------------
190      !!                    ***  ROUTINE sbc  ***
191      !!             
192      !! ** Purpose :   provide at each time-step the ocean surface boundary
193      !!                condition (momentum, heat and freshwater fluxes)
194      !!
195      !! ** Method  :   blah blah  to be written ?????????
196      !!                CAUTION : never mask the surface stress field (tke sbc)
197      !!
198      !! ** Action  : - set the ocean surface boundary condition at before and now
199      !!                time step, i.e. 
200      !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, emps_b, qrp_b, erp_b
201      !!                utau  , vtau  , qns  , qsr  , emp  , emps  , qrp  , erp
202      !!              - updte the ice fraction : fr_i
203      !!----------------------------------------------------------------------
204      INTEGER, INTENT(in) ::   kt       ! ocean time step
205      !!---------------------------------------------------------------------
206
207      !                                            ! ---------------------------------------- !
208      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
209         !                                         ! ---------------------------------------- !
210         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields
211         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
212         qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine)
213         ! The 3D heat content due to qsr forcing is treated in traqsr
214         ! qsr_b (:,:) = qsr (:,:)
215         emp_b (:,:) = emp (:,:)
216         emps_b(:,:) = emps(:,:)
217      ENDIF
218      !                                            ! ---------------------------------------- !
219      !                                            !        forcing field computation         !
220      !                                            ! ---------------------------------------- !
221
222      CALL iom_setkt( kt + nn_fsbc - 1 )                 ! in sbc, iom_put is called every nn_fsbc time step
223      !
224      IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc
225                                                         ! (caution called before sbc_ssm)
226      !
227      CALL sbc_ssm( kt )                                 ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
228      !                                                  ! averaged over nf_sbc time-step
229
230                                                   !==  sbc formulation  ==!
231                                                           
232      SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
233      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, emps)
234      CASE(  0 )   ;   CALL sbc_gyre    ( kt )                    ! analytical formulation : GYRE configuration
235      CASE(  1 )   ;   CALL sbc_ana     ( kt )                    ! analytical formulation : uniform sbc
236      CASE(  2 )   ;   CALL sbc_flx     ( kt )                    ! flux formulation
237      CASE(  3 )   ;   CALL sbc_blk_clio( kt )                    ! bulk formulation : CLIO for the ocean
238      CASE(  4 )   ;   CALL sbc_blk_core( kt )                    ! bulk formulation : CORE for the ocean
239      CASE(  5 )   ;   CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! coupled formulation
240      CASE( -1 )                               
241                       CALL sbc_ana     ( kt )                    ! ESOPA, test ALL the formulations
242                       CALL sbc_gyre    ( kt )                    !
243                       CALL sbc_flx     ( kt )                    !
244                       CALL sbc_blk_clio( kt )                    !
245                       CALL sbc_blk_core( kt )                    !
246                       CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   !
247      END SELECT
248
249      !                                            !==  Misc. Options  ==!
250     
251      SELECT CASE( nn_ice )                                     ! Update heat and freshwater fluxes over sea-ice areas
252      CASE(  1 )   ;       CALL sbc_ice_if   ( kt )                  ! Ice-cover climatology ("Ice-if" model)
253         !                                                     
254      CASE(  2 )   ;       CALL sbc_ice_lim_2( kt, nsbc )            ! LIM-2 ice model
255         IF( lk_obc )      CALL obc_ice_lim_2( kt )                  ! OBC boundary condition
256         !                                                     
257      CASE(  3 )   ;       CALL sbc_ice_lim  ( kt, nsbc )            ! LIM-3 ice model
258      END SELECT                                             
259
260      IF( ln_rnf       )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes
261 
262      IF( ln_ssr       )   CALL sbc_ssr( kt )                   ! add SST/SSS damping term
263
264      IF( nn_fwb  /= 0 )   CALL sbc_fwb( kt, nn_fwb, nn_fsbc )  ! control the freshwater budget
265
266      IF( nclosea == 1 )   CALL sbc_clo( kt )                   ! treatment of closed sea in the model domain
267      !                                                         ! (update freshwater fluxes)
268!RBbug do not understand why see ticket 667
269      CALL lbc_lnk( emp, 'T', 1. )
270      !
271      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
272         !                                             ! ---------------------------------------- !
273         IF( ln_rstart .AND.    &                               !* Restart: read in restart file
274            & iom_varid( numror, 'utau_b', ldstop = .FALSE. ) > 0 ) THEN
275            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file'
276            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b )   ! before i-stress  (U-point)
277            CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b )   ! before j-stress  (V-point)
278            CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b  )   ! before non solar heat flux (T-point)
279            ! The 3D heat content due to qsr forcing is treated in traqsr
280            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  )   ! before     solar heat flux (T-point)
281            CALL iom_get( numror, jpdom_autoglo, 'emp_b' , emp_b  )   ! before     freshwater flux (T-point)
282            CALL iom_get( numror, jpdom_autoglo, 'emps_b', emps_b )   ! before C/D freshwater flux (T-point)
283         ELSE                                                   !* no restart: set from nit000 values
284            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
285            utau_b(:,:) = utau(:,:) 
286            vtau_b(:,:) = vtau(:,:)
287            qns_b (:,:) = qns (:,:)
288            ! qsr_b (:,:) = qsr (:,:)
289            emp_b (:,:) = emp (:,:)
290            emps_b(:,:) = emps(:,:)
291         ENDIF
292      ENDIF
293      !                                                ! ---------------------------------------- !
294      IF( lrst_oce ) THEN                              !      Write in the ocean restart file     !
295         !                                             ! ---------------------------------------- !
296         IF(lwp) WRITE(numout,*)
297         IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   &
298            &                    'at it= ', kt,' date= ', ndastp
299         IF(lwp) WRITE(numout,*) '~~~~'
300         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau )
301         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau )
302         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  )
303         ! The 3D heat content due to qsr forcing is treated in traqsr
304         ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
305         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  )
306         CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emps )
307      ENDIF
308
309      !                                                ! ---------------------------------------- !
310      !                                                !        Outputs and control print         !
311      !                                                ! ---------------------------------------- !
312      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
313         CALL iom_put( "empmr" , emp  - rnf )                   ! upward water flux
314         CALL iom_put( "empsmr", emps - rnf )                   ! c/d water flux
315         CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux
316         CALL iom_put( "qns"   , qns        )                   ! solar heat flux
317         CALL iom_put( "qsr"   ,       qsr  )                   ! solar heat flux
318         IF( nn_ice > 0 )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction
319      ENDIF
320      !
321      CALL iom_setkt( kt )           ! iom_put outside of sbc is called at every time step
322      !
323      CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at
324      CALL iom_put( "vtau", vtau )   ! j-wind stress    each time step in sea-ice)
325      CALL iom_put( "taum", taum )   ! wind stress module
326      CALL iom_put( "wspd", wndm )   ! wind speed  module
327      !
328      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
329         CALL prt_ctl(tab2d_1=fr_i      , clinfo1=' fr_i     - : ', mask1=tmask, ovlap=1 )
330         CALL prt_ctl(tab2d_1=(emp-rnf) , clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 )
331         CALL prt_ctl(tab2d_1=(emps-rnf), clinfo1=' emps-rnf - : ', mask1=tmask, ovlap=1 )
332         CALL prt_ctl(tab2d_1=qns       , clinfo1=' qns      - : ', mask1=tmask, ovlap=1 )
333         CALL prt_ctl(tab2d_1=qsr       , clinfo1=' qsr      - : ', mask1=tmask, ovlap=1 )
334         CALL prt_ctl(tab3d_1=tmask     , clinfo1=' tmask    - : ', mask1=tmask, ovlap=1, kdim=jpk )
335         CALL prt_ctl(tab3d_1=tn        , clinfo1=' sst      - : ', mask1=tmask, ovlap=1, kdim=1   )
336         CALL prt_ctl(tab3d_1=sn        , clinfo1=' sss      - : ', mask1=tmask, ovlap=1, kdim=1   )
337         CALL prt_ctl(tab2d_1=utau      , clinfo1=' utau     - : ', mask1=umask,                      &
338            &         tab2d_2=vtau      , clinfo2=' vtau     - : ', mask2=vmask, ovlap=1 )
339      ENDIF
340      !
341   END SUBROUTINE sbc
342
343   !!======================================================================
344END MODULE sbcmod
Note: See TracBrowser for help on using the repository browser.