source: NEMO/branches/UKMO/NEMO_4.0_surge/cfgs/AMM15_SURGE/MY_SRC/usrdef_sbc.F90 @ 11180

Last change on this file since 11180 was 11180, checked in by clne, 2 years ago

Initial commit of code for 2d (surge) work in NEMO4.
This is aiming to replicate the 3.6 version in branches/UKMO/dev_r5518_Surge_Modelling

File size: 14.0 KB
Line 
1MODULE usrdef_sbc
2   !!======================================================================
3   !!                       ***  MODULE usrdef_sbc  ***
4   !!
5   !!                  ===  AMM7_SURGE configuration  ===
6   !!
7   !! User defined :   surface forcing of a user configuration
8   !!======================================================================
9   !! History :  4.0   ! 2016-03  (S. Flavoni, G. Madec)  user defined interface
10   !!            4.0   ! 2017-12  (C. O'Neill)  add necessary options for surge work - either no fluxes
11   !!                                           (for tide-only run) or wind and pressure only
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   usr_def_sbc    : user defined surface bounday conditions in LOCK_EXCHANGE case
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers
18   USE dom_oce         ! ocean space and time domain
19   USE sbc_oce         ! Surface boundary condition: ocean fields
20   USE sbc_ice         ! Surface boundary condition: ocean fields
21   USE fldread         ! read input fields
22   USE phycst          ! physical constants
23   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
24   !
25   USE in_out_manager  ! I/O manager
26   USE iom
27   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
28   USE lib_mpp         ! distribued memory computing library
29   !USE wrk_nemo       ! work arrays
30   USE timing         ! Timing
31   USE prtctl         ! Print control
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   usrdef_sbc_oce    ! routine called in sbcmod module
37   PUBLIC   usrdef_sbc_ice_tau  ! routine called by sbcice_lim.F90 for ice dynamics
38   PUBLIC   usrdef_sbc_ice_flx  ! routine called by sbcice_lim.F90 for ice thermo
39   PUBLIC   surge_oce    ! routine called by usrdef_sbc_oce (if required)
40   
41   
42   INTEGER , PARAMETER ::   jpfld   = 2           ! maximum number of files to read
43   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point
44   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point
45
46   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
47
48   REAL(wp), PARAMETER ::   rhoa =    1.22        ! air density
49
50   !                                  !!* Namelist namsbc_usr
51   REAL(wp) ::   rn_vfac     ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
52   REAL(wp) ::   rn_charn_const 
53   LOGICAL  ::   ln_use_sbc  ! Surface fluxes on or not
54
55
56   !! * Substitutions
57#  include "vectopt_loop_substitute.h90"
58   !!----------------------------------------------------------------------
59   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
60   !! $Id$
61   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
62   !!----------------------------------------------------------------------
63CONTAINS
64
65   SUBROUTINE usrdef_sbc_oce( kt )
66      !!---------------------------------------------------------------------
67      !!                    ***  ROUTINE usr_def_sbc  ***
68      !!             
69      !! ** Purpose :   provide at each time-step the surface boundary
70      !!              condition, i.e. the momentum, heat and freshwater fluxes.
71      !!
72      !! ** Method  :   all 0 fields, for AMM7_SURGE case
73      !!                CAUTION : never mask the surface stress field !
74      !!
75      !! ** Action  : - if tide-only case - set to ZERO all the ocean surface boundary condition, i.e.   
76      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
77      !!              - if tide+surge case - read in wind and air pressure      !!
78      !!----------------------------------------------------------------------
79      INTEGER, INTENT(in) ::   kt   ! ocean time step
80
81      INTEGER  ::   ierror   ! return error code
82      INTEGER  ::   ifpr     ! dummy loop indice
83      INTEGER  ::   ios      ! Local integer output status for namelist read
84      !
85      CHARACTER(len=100) ::  cn_dir   !   Root directory for location of flux files
86      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read
87      TYPE(FLD_N) ::   sn_wndi, sn_wndj                        ! informations about the fields to be read
88
89      NAMELIST/namsbc_usr/ ln_use_sbc, cn_dir , rn_vfac,  &
90         &                   sn_wndi, sn_wndj, rn_charn_const
91      !!---------------------------------------------------------------------
92      !
93      IF( kt == nit000 ) THEN
94         
95         
96         REWIND( numnam_cfg )              ! Namelist namsbc_usr in configuration namelist
97         READ  ( numnam_cfg, namsbc_usr, IOSTAT = ios, ERR = 902 )
98902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_surge in configuration namelist', lwp )
99
100         IF(lwm) WRITE( numond, namsbc_usr )
101
102         IF(ln_use_sbc) THEN
103             IF(lwp) WRITE(numout,*)' usr_sbc : AMM7_SURGE tide + surge case: surface wind and pressure (assuming ln_dyn_apr=T) applied'
104             IF(lwp) WRITE(numout,*)' ~~~~~~~~~~~ '
105
106             !                                         ! store namelist information in an array
107             slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
108             !
109             !
110             ALLOCATE( sf(jpfld), STAT=ierror )         ! set sf structure
111             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_surge: unable to allocate sf structure' )
112             DO ifpr= 1, jpfld
113                ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
114                IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
115             END DO
116             !                                         ! fill sf with slf_i and control print
117             CALL fld_fill( sf, slf_i, cn_dir, 'sbc_surge', 'flux formulation for ocean surface boundary condition', 'namsbc_surge' )
118
119         ELSE
120             IF(lwp) WRITE(numout,*)' usr_sbc : AMM7_SURGE tide only case: NO surface forcing'
121             IF(lwp) WRITE(numout,*)' ~~~~~~~~~~~   utau = vtau = taum = wndm = qns = qsr = emp = sfx = 0'
122
123             utau(:,:) = 0._wp
124             vtau(:,:) = 0._wp
125             taum(:,:) = 0._wp
126             wndm(:,:) = 0._wp
127             !
128             emp (:,:) = 0._wp
129             sfx (:,:) = 0._wp
130             qns (:,:) = 0._wp
131             qsr (:,:) = 0._wp
132             !         
133             uwnd(:,:) = 0._wp
134             vwnd(:,:) = 0._wp
135         ENDIF
136
137      ENDIF
138      !
139      IF(ln_use_sbc) THEN
140          CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
141
142          !                                            ! compute the surface ocean fluxes using CORE bulk formulea
143          IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL surge_oce( kt, sf, ssu_m, ssv_m, rn_charn_const )
144         
145      ENDIF
146   END SUBROUTINE usrdef_sbc_oce
147
148
149   
150   SUBROUTINE surge_oce( kt, sf, pu, pv, rn_charn_const )
151      !!---------------------------------------------------------------------
152      !!                     ***  ROUTINE surge_oce  ***
153      !!
154      !! ** Purpose :   provide the momentum fluxes at
155      !!      the ocean surface at each time step
156      !!
157      !! ** Method  :   Charnock formulea for the ocean using atmospheric
158      !!      fields read in sbc_read
159      !!
160      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
161      !!              - vtau    : j-component of the stress at V-point  (N/m2)
162      !!              - taum    : Wind stress module at T-point         (N/m2)
163      !!              - wndm    : Wind speed module at T-point          (m/s)
164      !!
165      !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC
166      !!---------------------------------------------------------------------
167      INTEGER  , INTENT(in   )                 ::   kt    ! time step index
168      TYPE(fld), INTENT(inout), DIMENSION(:)   ::   sf    ! input data
169      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pu    ! surface current at U-point (i-component) [m/s]
170      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pv    ! surface current at V-point (j-component) [m/s]
171      REAL(wp) , INTENT(in)                    ::   rn_charn_const! local variable
172      !
173      INTEGER  ::   ji, jj               ! dummy loop indices
174      REAL(wp) ::   zztmp                ! local variable
175      REAL(wp) ::   z_z0, z_Cd1          ! local variable
176      REAL(wp) ::   zi                   ! local variable
177      REAL(wp), ALLOCATABLE, DIMENSION(:,:)::   zwnd_i, zwnd_j    ! wind speed components at T-point
178      REAL(wp), ALLOCATABLE, DIMENSION(:,:)::   Cd                ! transfer coefficient for momentum      (tau)
179      !!---------------------------------------------------------------------
180      !
181      IF( ln_timing )  CALL timing_start('surge_oce')
182      !
183      ALLOCATE( zwnd_i(jpi,jpj) )
184      ALLOCATE( zwnd_j(jpi,jpj) )
185      ALLOCATE( Cd(jpi,jpj) )
186      !
187      ! ----------------------------------------------------------------------------- !
188      !      0   Wind components and module at T-point relative to the moving ocean   !
189      ! ----------------------------------------------------------------------------- !
190
191      ! ... components ( U10m - U_oce ) at T-point (unmasked)
192      zwnd_i(:,:) = 0.e0 
193      zwnd_j(:,:) = 0.e0
194      DO jj = 2, jpjm1
195         DO ji = fs_2, fs_jpim1   ! vect. opt.
196            uwnd(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
197            vwnd(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
198         END DO
199      END DO
200      zwnd_i(:,:) = uwnd(:,:)
201      zwnd_j(:,:) = vwnd(:,:)
202
203      CALL lbc_lnk( 'surge_oce', zwnd_i(:,:) , 'T', -1. )
204      CALL lbc_lnk( 'surge_oce', zwnd_j(:,:) , 'T', -1. )
205      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
206      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
207         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
208
209      ! ----------------------------------------------------------------------------- !
210      !      I   Radiative FLUXES                                                     !
211      ! ----------------------------------------------------------------------------- !
212     
213      qsr(:,:)=0._wp
214
215      ! ----------------------------------------------------------------------------- !
216      !     II    Turbulent FLUXES                                                    !
217      ! ----------------------------------------------------------------------------- !
218      Cd(:,:)=0.0001_wp
219      DO jj = 1,jpj
220         DO ji = 1,jpi
221            z_Cd1=0._wp
222            zi=1
223            !Iterate
224            DO WHILE((abs(Cd(ji,jj)-z_Cd1))>1E-6)
225               z_Cd1=Cd(ji,jj)
226               z_z0=rn_charn_const*z_Cd1*wndm(ji,jj)**2/grav
227               Cd(ji,jj)=(0.41_wp/log(10._wp/z_z0))**2
228               zi=zi+1
229            ENDDO
230         ENDDO
231      ENDDO
232
233      ! ... tau module, i and j component
234      DO jj = 1, jpj
235         DO ji = 1, jpi
236            zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj)
237            taum  (ji,jj) = zztmp * wndm  (ji,jj)
238            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
239            zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
240         END DO
241      END DO
242
243      CALL iom_put( "taum_oce", taum )   ! output wind stress module
244      CALL iom_put( "uwnd", uwnd )   ! output wind stress module
245      CALL iom_put( "vwnd", vwnd )   ! output wind stress module
246
247      ! ... utau, vtau at U- and V_points, resp.
248      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
249      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
250      DO jj = 1, jpjm1
251         DO ji = 1, fs_jpim1
252            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) &
253               &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
254            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) &
255               &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
256         END DO
257      END DO
258      CALL lbc_lnk( 'surge_oce', utau(:,:), 'U', -1. )
259      CALL lbc_lnk( 'surge_oce', vtau(:,:), 'V', -1. )
260
261   
262      IF(ln_ctl) THEN
263         CALL prt_ctl( tab2d_1=utau  , clinfo1=' surge_oce: utau   : ', mask1=umask,   &
264            &          tab2d_2=vtau  , clinfo2=           ' vtau : '  , mask2=vmask )
265         CALL prt_ctl( tab2d_1=wndm  , clinfo1=' surge_oce: wndm   : ')
266      ENDIF
267       
268      ! ----------------------------------------------------------------------------- !
269      !     III    Total FLUXES                                                       !
270      ! ----------------------------------------------------------------------------- !
271      !
272      emp (:,:) = 0._wp
273      qns(:,:)  = 0._wp
274      sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
275      !
276!       IF ( nn_ice == 0 ) THEN
277!          CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean
278!          CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean
279!          CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean
280!       ENDIF
281      !
282      IF(ln_ctl) THEN
283         CALL prt_ctl(tab2d_1=utau , clinfo1=' surge_oce: utau   : ', mask1=umask,   &
284            &         tab2d_2=vtau , clinfo2=           ' vtau  : ' , mask2=vmask )
285      ENDIF
286      !
287!       DEALLOCATE( zwnd_i(jpi,jpj) )
288!       DEALLOCATE( zwnd_j(jpi,jpj) )
289!       DEALLOCATE( Cd(jpi,jpj) )
290      !
291      IF( ln_timing )  CALL timing_stop('surge_oce')
292
293      !
294   END SUBROUTINE surge_oce
295
296
297   SUBROUTINE usrdef_sbc_ice_tau( kt )
298      INTEGER, INTENT(in) ::   kt   ! ocean time step
299   END SUBROUTINE usrdef_sbc_ice_tau
300
301   SUBROUTINE usrdef_sbc_ice_flx( kt )
302      INTEGER, INTENT(in) ::   kt   ! ocean time step
303   END SUBROUTINE usrdef_sbc_ice_flx
304
305   !!======================================================================
306END MODULE usrdef_sbc
Note: See TracBrowser for help on using the repository browser.