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 @ 7168

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

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