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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90 @ 4416

Last change on this file since 4416 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

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