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.
sbcwave.F90 in branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90 @ 7193

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

Update latest changes to the INGV branch, fix small bug (read ln_zdfqiao from the namelist)
File size: 16.9 KB
Line 
1MODULE sbcwave
2   !!======================================================================
3   !!                       ***  MODULE  sbcwave  ***
4   !! Wave module
5   !!======================================================================
6   !! History :  3.3  !   2011-09  (Adani M)  Original code: Drag Coefficient
7   !!         :  3.4  !   2012-10  (Adani M)                 Stokes Drift
8   !!            3.6  !   2014-09  (Clementi E, Oddo P)New Stokes Drift Computation
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   sbc_wave      : wave data from wave model in netcdf files
13   !!----------------------------------------------------------------------
14   USE oce            !
15   USE sbc_oce       ! Surface boundary condition: ocean fields
16   USE bdy_oce        !
17   USE domvvl         !
18   !
19   USE iom            ! I/O manager library
20   USE in_out_manager ! I/O manager
21   USE lib_mpp        ! distribued memory computing library
22   USE fldread       ! read input fields
23   USE wrk_nemo       !
24   USE phycst         ! physical constants
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   sbc_stokes, sbc_qiao  ! routines called in sbccpl
30   PUBLIC   sbc_wave    ! routine called in sbcmod
31   
32   ! Variables checking if the wave parameters are coupled (if not, they are read from file)
33   LOGICAL, PUBLIC     ::   cpl_hsig=.FALSE.
34   LOGICAL, PUBLIC     ::   cpl_phioc=.FALSE.
35   LOGICAL, PUBLIC     ::   cpl_sdrftx=.FALSE.
36   LOGICAL, PUBLIC     ::   cpl_sdrfty=.FALSE.
37   LOGICAL, PUBLIC     ::   cpl_wper=.FALSE.
38   LOGICAL, PUBLIC     ::   cpl_wnum=.FALSE.
39   LOGICAL, PUBLIC     ::   cpl_wstrf=.FALSE.
40   LOGICAL, PUBLIC     ::   cpl_wdrag=.FALSE.
41
42   INTEGER             ::   jpfld                ! number of files to read for stokes drift
43   INTEGER             ::   jp_usd               ! index of stokes drift  (i-component) (m/s)    at T-point
44   INTEGER             ::   jp_vsd               ! index of stokes drift  (j-component) (m/s)    at T-point
45   INTEGER             ::   jp_swh               ! index of significant wave hight      (m)      at T-point
46   INTEGER             ::   jp_wmp               ! index of mean wave period            (s)      at T-point
47
48   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient
49   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift
50   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao
51   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean
52   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: cdn_wave 
53   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: swh,wmp, wnum
54   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tauoc_wave
55   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tsd2d
56   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: usd2d, vsd2d
57   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     :: usd3d, vsd3d, wsd3d 
58   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     :: usd3dt, vsd3dt
59   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: zusd2dt, zvsd2dt
60
61   !! * Substitutions
62#  include "domzgr_substitute.h90"
63#  include "vectopt_loop_substitute.h90"
64   !!----------------------------------------------------------------------
65   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
66   !! $Id$
67   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
68   !!----------------------------------------------------------------------
69CONTAINS
70
71   SUBROUTINE sbc_stokes( )
72      !!---------------------------------------------------------------------
73      !!                     ***  ROUTINE sbc_stokes  ***
74      !!
75      !! ** Purpose :   compute the 3d Stokes Drift according to Breivik et al.,
76      !!                2014 (DOI: 10.1175/JPO-D-14-0020.1)
77      !!
78      !! ** Method  : - Calculate Stokes transport speed
79      !!              - Calculate horizontal divergence
80      !!              - Integrate the horizontal divergenze from the bottom
81      !! ** action 
82      !!---------------------------------------------------------------------
83      INTEGER                ::   jj,ji,jk 
84      REAL(wp)                       ::  ztransp, zsp0, zk, zus, zvs
85      REAL(wp), DIMENSION(jpi,jpj)   ::  zfac 
86      REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv   ! 3D workspace
87
88
89      CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv )
90      DO jk = 1, jpk
91         DO jj = 1, jpj
92            DO ji = 1, jpi
93            ! On T grid
94            ! Stokes transport speed estimated from Hs and Tmean
95            ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp))
96            ! Stokes surface speed
97            zsp0 = SQRT( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2)
98            ! Wavenumber scale
99            zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp)
100            ! Depth attenuation
101            zfac(ji,jj) = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk))
102            END DO
103         END DO
104         !
105         DO jj = 1, jpj
106            DO ji = 1, jpi
107               usd3dt(ji,jj,jk) = zfac(ji,jj) * zusd2dt(ji,jj) * tmask(ji,jj,jk) 
108               vsd3dt(ji,jj,jk) = zfac(ji,jj) * zvsd2dt (ji,jj) * tmask(ji,jj,jk) 
109            END DO
110         END DO
111      END DO 
112      ! Into the U and V Grid
113      DO jk = 1, jpkm1 
114         DO jj = 1, jpjm1
115            DO ji = 1, jpim1
116               usd3d(ji,jj,jk) =   0.5 *  umask(ji,jj,jk)  *        &
117                               &  (usd3dt(ji,jj,jk) + usd3dt(ji+1,jj,jk))
118               vsd3d(ji,jj,jk) =  0.5 *  vmask(ji,jj,jk)  *           &
119                               &  (vsd3dt(ji,jj,jk) + vsd3dt(ji,jj+1,jk))
120            END DO
121         END DO
122      END DO
123      !
124      CALL lbc_lnk( usd3d(:,:,:), 'U', -1. )
125      CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. )
126      !
127      DO jk = 1, jpkm1                       ! e3t * Horizontal divergence
128         DO jj = 2, jpjm1
129            DO ji = fs_2, fs_jpim1   ! vector opt.
130               ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * fse3u_n(ji  ,jj,jk) * usd3d(ji  ,jj,jk)     &
131                  &                 - e2u(ji-1,jj) * fse3u_n(ji-1,jj,jk) * usd3d(ji-1,jj,jk)     &
132                  &                 + e1v(ji,jj  ) * fse3v_n(ji,jj  ,jk) * vsd3d(ji,jj  ,jk)     &
133                  &                 - e1v(ji,jj-1) * fse3v_n(ji,jj-1,jk) * vsd3d(ji,jj-1,jk)   ) * r1_e1e2t(ji,jj)
134            END DO
135         END DO
136         IF( .NOT. AGRIF_Root() ) THEN
137            IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   : ,jk) = 0._wp      ! east
138            IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   : ,jk) = 0._wp      ! west
139            IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :  ,nlcj-1,jk) = 0._wp      ! north
140            IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2  ,jk) = 0._wp      ! south
141         ENDIF
142      END DO
143      CALL lbc_lnk( ze3hdiv, 'T', 1. )
144      !
145      DO jk = jpkm1, 1, -1                   !* integrate from the bottom the e3t * hor. divergence
146         wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - ze3hdiv(:,:,jk)
147      END DO
148#if defined key_bdy
149      IF( lk_bdy ) THEN
150         DO jk = 1, jpkm1
151            wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:)
152         END DO
153      ENDIF
154#endif
155      CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv )
156      !
157   END SUBROUTINE sbc_stokes
158
159   SUBROUTINE sbc_qiao( )
160      !!---------------------------------------------------------------------
161      !!                     ***  ROUTINE sbc_qiao  ***
162      !!
163      !! ** Purpose :   Qiao formulation for wave enhanced turbulence
164      !!                2010 (DOI: 10.1007/s10236-010-0326)
165      !!
166      !! ** Method  : -
167      !! ** action 
168      !!---------------------------------------------------------------------
169      INTEGER                ::   jj,ji
170
171      ! Calculate the module of the stokes drift on T grid
172      !-------------------------------------------------
173      DO jj = 1, jpj
174         DO ji = 1, jpi
175            tsd2d(ji,jj) = SQRT( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2)
176         END DO
177      END DO
178      !
179   END SUBROUTINE sbc_qiao
180
181   SUBROUTINE sbc_wave( kt )
182      !!---------------------------------------------------------------------
183      !!                     ***  ROUTINE sbc_wave  ***
184      !!
185      !! ** Purpose :   read wave parameters from wave model  in netcdf files.
186      !!
187      !! ** Method  : - Read namelist namsbc_wave
188      !!              - Read Cd_n10 fields in netcdf files
189      !!              - Read stokes drift 2d in netcdf files
190      !!              - Read wave number in netcdf files
191      !!              - Compute 3d stokes drift using Breivik et al.,2014
192      !!                formulation
193      !! ** action 
194      !!---------------------------------------------------------------------
195      USE zdf_oce    ,  ONLY : ln_zdfqiao
196
197      IMPLICIT NONE
198
199      INTEGER, INTENT( in  ) ::   kt       ! ocean time step
200      !
201      INTEGER                ::   ierror   ! return error code
202      INTEGER                ::   ifpr
203      INTEGER                ::   ios      ! Local integer output status for namelist read
204      !
205      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files
206      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i     ! array of namelist informations on the fields to read
207      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  &
208                             &   sn_swh, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read
209      !!
210      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_swh, sn_wmp, sn_wnum, sn_tauoc
211      !!---------------------------------------------------------------------
212      !
213      !                                         ! -------------------- !
214      IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
215         !                                      ! -------------------- !
216         REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model
217         READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901)
218901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp )
219
220         REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model
221         READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 )
222902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp )
223         IF(lwm) WRITE ( numond, namsbc_wave )
224         !
225         IF ( ln_cdgw ) THEN
226            IF ( .NOT. cpl_wdrag ) THEN
227               ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
228               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
229               !
230                                      ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   )
231               IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) )
232               CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
233            ENDIF
234            ALLOCATE( cdn_wave(jpi,jpj) )
235            cdn_wave(:,:) = 0.0
236         ENDIF
237
238         IF ( ln_tauoc ) THEN
239            IF ( .NOT. cpl_wstrf ) THEN
240               ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc
241               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
242               !
243                                       ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   )
244               IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) )
245               CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
246            ENDIF
247            ALLOCATE( tauoc_wave(jpi,jpj) )
248            tauoc_wave(:,:) = 0.0
249         ENDIF
250
251         IF ( ln_sdw ) THEN
252            ! Find out how many fields have to be read from file if not coupled
253            jpfld=0
254            jp_usd=0; jp_vsd=0; jp_swh=0; jp_wmp=0
255            IF( .NOT. cpl_sdrftx ) THEN
256               jpfld=jpfld+1
257               jp_usd=jpfld
258            ENDIF
259            IF( .NOT. cpl_sdrfty ) THEN
260               jpfld=jpfld+1
261               jp_vsd=jpfld
262            ENDIF
263            IF( .NOT. cpl_hsig ) THEN
264               jpfld=jpfld+1
265               jp_swh=jpfld
266            ENDIF
267            IF( .NOT. cpl_wper ) THEN
268               jpfld=jpfld+1
269               jp_wmp=jpfld
270            ENDIF
271
272            ! Read from file only the non-coupled fields
273            IF( jpfld > 0 ) THEN
274               ALLOCATE( slf_i(jpfld) )
275               IF( jp_usd > 0 ) slf_i(jp_usd) = sn_usd
276               IF( jp_vsd > 0 ) slf_i(jp_vsd) = sn_vsd
277               IF( jp_swh > 0 ) slf_i(jp_swh) = sn_swh
278               IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp
279               ALLOCATE( sf_sd(jpfld), STAT=ierror )           !* allocate and fill sf_sd with stokes drift
280               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
281               !
282               DO ifpr= 1, jpfld
283                  ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) )
284                  IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) )
285               END DO
286
287               CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
288            ENDIF
289            ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj) )
290            ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) )
291            ALLOCATE( usd3dt(jpi,jpj,jpk),vsd3dt(jpi,jpj,jpk) )
292            ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) )
293            ALLOCATE( zusd2dt(jpi,jpj), zvsd2dt(jpi,jpj) )
294            usd3d(:,:,:) = 0._wp   ;   usd2d(:,:) = 0._wp   ;   
295            vsd3d(:,:,:) = 0._wp   ;   vsd2d(:,:) = 0._wp   ; 
296            wsd3d(:,:,:) = 0._wp   ;
297            usd3dt(:,:,:) = 0._wp  ;   vsd3dt(:,:,:) = 0._wp   ;
298            swh  (:,:)   = 0._wp   ;    wmp (:,:) = 0._wp ;
299            IF ( ln_zdfqiao ) THEN     !==  Vertical mixing enhancement using Qiao,2010  ==!
300               IF ( .NOT. cpl_wnum ) THEN
301                  ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
302                  IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable toallocate sf_wave structure' )
303                                         ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   )
304                  IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) )
305                  CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
306               ENDIF
307               ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) )
308               wnum(:,:) = 0._wp ; tsd2d(:,:) = 0._wp
309            ENDIF
310         ENDIF
311      ENDIF
312      !
313      IF ( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN              !==  Neutral drag coefficient  ==!
314         CALL fld_read( kt, nn_fsbc, sf_cd )      ! read from external forcing
315         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
316      ENDIF
317
318      IF ( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN             !==  Wave induced stress  ==!
319         CALL fld_read( kt, nn_fsbc, sf_tauoc )      !* read wave norm stress from external forcing
320         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1)
321      ENDIF
322
323      IF ( ln_sdw )  THEN                         !==  Computation of the 3d Stokes Drift  ==!
324         !
325         ! Read from file only if the field is not coupled
326         IF( jpfld > 0 ) THEN
327            CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing
328            IF( jp_swh > 0 ) swh(:,:)     = sf_sd(jp_swh)%fnow(:,:,1)   ! significant wave height
329            IF( jp_wmp > 0 ) wmp(:,:)     = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period
330            IF( jp_usd > 0 ) zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point
331            IF( jp_vsd > 0 ) zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point
332         ENDIF
333         !
334         ! Read also wave number if needed, so that it is available in coupling routines
335         IF ( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN
336            CALL fld_read( kt, nn_fsbc, sf_wn )      !* read wave parameters from external forcing
337            wnum(:,:) = sf_wn(1)%fnow(:,:,1)
338         ENDIF
339           
340         !==  Computation of the 3d Stokes Drift according to Breivik et al.,2014
341         !(DOI: 10.1175/JPO-D-14-0020.1)==!
342         !
343         ! Calculate only if no necessary fields are coupled, if not calculate later after coupling
344         IF( jpfld == 4 ) THEN
345            CALL sbc_stokes()
346            IF ( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN
347               CALL sbc_qiao()
348            ENDIF
349         ENDIF
350      ENDIF
351      !
352   END SUBROUTINE sbc_wave
353     
354   !!======================================================================
355END MODULE sbcwave
Note: See TracBrowser for help on using the repository browser.