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.
sbcsurge.F90 in branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcsurge.F90 @ 6592

Last change on this file since 6592 was 6592, checked in by rfurner, 8 years ago

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