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

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90 @ 9115

Last change on this file since 9115 was 9115, checked in by acc, 6 years ago

Branch 2017/dev_merge_2017. Tidier version of sbcwave.F90 with associated changes in modules. References and schemes for Stokes drift parameterisations have been updated. Choices now are the original Breivik 2014 scheme and the latest Li et al 2017 scheme (which is based on the Breivik et al 2016 Phillips spectrum but with a depth averaged profile). This has been compiled but not yet tested. Also removed trailing spaces in lbc_lnk_multi_generic.h90
  • Property svn:keywords set to Id
File size: 25.2 KB
Line 
1MODULE sbcwave
2   !!======================================================================
3   !!                       ***  MODULE  sbcwave  ***
4   !! Wave module
5   !!======================================================================
6   !! History :  3.3  !  2011-09  (M. Adani)  Original code: Drag Coefficient
7   !!         :  3.4  !  2012-10  (M. Adani)  Stokes Drift
8   !!            3.6  !  2014-09  (E. Clementi,P. Oddo) New Stokes Drift Computation
9   !!             -   !  2016-12  (G. Madec, E. Clementi) update Stoke drift computation
10   !!                                                    + add sbc_wave_ini routine
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   sbc_stokes    : calculate 3D Stokes-drift velocities
15   !!   sbc_wave      : wave data from wave model in netcdf files
16   !!   sbc_wave_init : initialisation fo surface waves
17   !!----------------------------------------------------------------------
18   USE phycst         ! physical constants
19   USE oce            ! ocean variables
20   USE sbc_oce        ! Surface boundary condition: ocean fields
21   USE zdf_oce,  ONLY : ln_zdfswm
22   USE bdy_oce        ! open boundary condition variables
23   USE domvvl         ! domain: variable volume layers
24   !
25   USE iom            ! I/O manager library
26   USE in_out_manager ! I/O manager
27   USE lib_mpp        ! distribued memory computing library
28   USE fldread        ! read input fields
29   USE wrk_nemo       !
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   sbc_stokes      ! routine called in sbccpl
35   PUBLIC   sbc_wstress     ! routine called in sbcmod
36   PUBLIC   sbc_wave        ! routine called in sbcmod
37   PUBLIC   sbc_wave_init   ! routine called in sbcmod
38   
39   ! Variables checking if the wave parameters are coupled (if not, they are read from file)
40   LOGICAL, PUBLIC ::   cpl_hsig   = .FALSE.
41   LOGICAL, PUBLIC ::   cpl_phioc  = .FALSE.
42   LOGICAL, PUBLIC ::   cpl_sdrftx = .FALSE.
43   LOGICAL, PUBLIC ::   cpl_sdrfty = .FALSE.
44   LOGICAL, PUBLIC ::   cpl_wper   = .FALSE.
45   LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE.
46   LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE.
47   LOGICAL, PUBLIC ::   cpl_tauwoc = .FALSE.
48   LOGICAL, PUBLIC ::   cpl_tauw   = .FALSE.
49   LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE.
50
51   INTEGER ::   jpfld    ! number of files to read for stokes drift
52   INTEGER ::   jp_usd   ! index of stokes drift  (i-component) (m/s)    at T-point
53   INTEGER ::   jp_vsd   ! index of stokes drift  (j-component) (m/s)    at T-point
54   INTEGER ::   jp_hsw   ! index of significant wave hight      (m)      at T-point
55   INTEGER ::   jp_wmp   ! index of mean wave period            (s)      at T-point
56   INTEGER ::   jp_wfr   ! index of wave peak frequency         (1/s)    at T-point
57
58   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_cd      ! structure of input fields (file informations, fields read) Drag Coefficient
59   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sd      ! structure of input fields (file informations, fields read) Stokes Drift
60   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_wn      ! structure of input fields (file informations, fields read) wave number for Qiao
61   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tauwoc  ! structure of input fields (file informations, fields read) normalized wave stress into the ocean
62   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tauw    ! structure of input fields (file informations, fields read) ocean stress components from wave model
63
64   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !:
65   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw, wmp, wnum      !:
66   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   wfreq               !:
67   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !: 
68   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauw_x, tauw_y      !: 
69   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !:
70   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   div_sd              !: barotropic stokes drift divergence
71   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   ut0sd, vt0sd        !: surface Stokes drift velocities at t-point
72   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd  , vsd  , wsd   !: Stokes drift velocities at u-, v- & w-points, resp.
73
74   !! * Substitutions
75#  include "vectopt_loop_substitute.h90"
76   !!----------------------------------------------------------------------
77   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
78   !! $Id$
79   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
80   !!----------------------------------------------------------------------
81CONTAINS
82
83   SUBROUTINE sbc_stokes( )
84      !!---------------------------------------------------------------------
85      !!                     ***  ROUTINE sbc_stokes  ***
86      !!
87      !! ** Purpose :   compute the 3d Stokes Drift according to Breivik et al.,
88      !!                2014 (DOI: 10.1175/JPO-D-14-0020.1)
89      !!
90      !! ** Method  : - Calculate Stokes transport speed
91      !!              - Calculate horizontal divergence
92      !!              - Integrate the horizontal divergenze from the bottom
93      !! ** action 
94      !!---------------------------------------------------------------------
95      INTEGER  ::   jj, ji, jk   ! dummy loop argument
96      INTEGER  ::   ik           ! local integer
97      REAL(wp) ::  ztransp, zfac, zsp0
98      REAL(wp) ::  zdepth, zsqrt_depth,  zexp_depth, z_two_thirds, zsqrtpi !sqrt of pi
99      REAL(wp) ::  zbot_u, zbot_v, zkb_u, zkb_v, zke3_u, zke3_v, zda_u, zda_v
100      REAL(wp) ::  zstokes_psi_u_bot, zstokes_psi_v_bot
101      REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v
102      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd     ! 2D workspace
103      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zstokes_psi_u_top, zstokes_psi_v_top ! 2D workspace
104      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ze3divh                              ! 3D workspace
105      !!---------------------------------------------------------------------
106      !
107      ALLOCATE( ze3divh(jpi,jpj,jpk) )
108      ALLOCATE( zk_t(jpi,jpj), zk_u(jpi,jpj), zk_v(jpi,jpj), zu0_sd(jpi,jpj), zv0_sd(jpi,jpj) )
109      !
110      ! select parameterization for the calculation of vertical Stokes drift
111      ! exp. wave number at t-point
112      IF( ll_st_bv_li ) THEN   ! (Eq. (19) in Breivik et al. (2014) )
113         zfac = 2.0_wp * rpi / 16.0_wp
114         DO jj = 1, jpj
115            DO ji = 1, jpi
116               ! Stokes drift velocity estimated from Hs and Tmean
117               ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp )
118               ! Stokes surface speed
119               tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj))
120               ! Wavenumber scale
121               zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp )
122            END DO
123         END DO
124         DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points
125            DO ji = 1, jpim1
126               zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
127               zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
128               !
129               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
130               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
131            END DO
132         END DO
133      ELSE IF( ll_st_peakfr ) THEN    ! peak wave number calculated from the peak frequency received by the wave model
134         DO jj = 1, jpjm1
135            DO ji = 1, jpim1
136               zk_u(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji+1,jj)*wfreq(ji+1,jj) ) / grav
137               zk_v(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji,jj+1)*wfreq(ji,jj+1) ) / grav
138               !
139               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
140               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
141            END DO
142         END DO
143      ENDIF
144      !
145      !                       !==  horizontal Stokes Drift 3D velocity  ==!
146      IF( ll_st_bv2014 ) THEN
147         DO jk = 1, jpkm1
148            DO jj = 2, jpjm1
149               DO ji = 2, jpim1
150                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )
151                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )
152                  !                         
153                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
154                  zkh_v = zk_v(ji,jj) * zdep_v
155                  !                                ! Depth attenuation
156                  zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u )
157                  zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v )
158                  !
159                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
160                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
161               END DO
162            END DO
163         END DO
164      ELSE IF( ll_st_li2017 .OR. ll_st_peakfr ) THEN
165         ALLOCATE( zstokes_psi_u_top(jpi,jpj), zstokes_psi_v_top(jpi,jpj) )
166         DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points
167            DO ji = 1, jpim1
168               zstokes_psi_u_top(ji,jj) = 0._wp
169               zstokes_psi_v_top(ji,jj) = 0._wp
170            END DO
171         END DO
172         DO jk = 1, jpkm1
173            DO jj = 2, jpjm1
174               DO ji = 2, jpim1
175                  zbot_u = 0.5_wp * ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji+1,jj,jk+1) )
176                  zbot_v = 0.5_wp * ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji,jj+1,jk+1) )
177                  zkb_u  = 2.0_wp * zk_u(ji,jj) * zbot_u                             ! 2k * bottom depth
178                  zkb_v  = 2.0_wp * zk_v(ji,jj) * zbot_v                             ! 2k * bottom depth
179                  !
180                  zke3_u = MAX(1.e-8_wp, 2.0_wp * zk_u(ji,jj) * e3u_n(ji,jj,jk))     ! 2k * thickness
181                  zke3_v = MAX(1.e-8_wp, 2.0_wp * zk_v(ji,jj) * e3v_n(ji,jj,jk))     ! 2k * thickness
182
183                  ! Depth attenuation .... do u component first..
184                  zdepth      = zkb_u
185                  zsqrt_depth = SQRT(zdepth)
186                  zexp_depth  = EXP(-zdepth)
187                  zstokes_psi_u_bot = 1.0_wp - zexp_depth  &
188                       &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
189                       &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )
190                  zda_u                    = ( zstokes_psi_u_bot - zstokes_psi_u_top(ji,jj) ) / zke3_u
191                  zstokes_psi_u_top(ji,jj) =   zstokes_psi_u_bot
192
193                  !         ... and then v component
194                  zdepth      =zkb_v
195                  zsqrt_depth = SQRT(zdepth)
196                  zexp_depth  = EXP(-zdepth)
197                  zstokes_psi_v_bot = 1.0_wp - zexp_depth  &
198                       &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
199                       &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )
200                  zda_v                    = ( zstokes_psi_v_bot - zstokes_psi_v_top(ji,jj) ) / zke3_v
201                  zstokes_psi_v_top(ji,jj) =   zstokes_psi_v_bot
202                  !
203                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
204                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
205               END DO
206            END DO
207         END DO
208         DEALLOCATE( zstokes_psi_u_top, zstokes_psi_v_top )
209      ENDIF
210
211      CALL lbc_lnk_multi( usd, 'U', -1., vsd, 'V', -1. )
212
213      !
214      !                       !==  vertical Stokes Drift 3D velocity  ==!
215      !
216      DO jk = 1, jpkm1               ! Horizontal e3*divergence
217         DO jj = 2, jpj
218            DO ji = fs_2, jpi
219               ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * usd(ji  ,jj,jk)    &
220                  &                 - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * usd(ji-1,jj,jk)    &
221                  &                 + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vsd(ji,jj  ,jk)    &
222                  &                 - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vsd(ji,jj-1,jk)  ) * r1_e1e2t(ji,jj)
223            END DO
224         END DO
225      END DO
226      !
227#if defined key_agrif
228      IF( .NOT. Agrif_Root() ) THEN
229         IF( nbondi == -1 .OR. nbondi == 2 )   ze3divh( 2:nbghostcells+1,:        ,:) = 0._wp      ! west
230         IF( nbondi ==  1 .OR. nbondi == 2 )   ze3divh( nlci-nbghostcells:nlci-1,:,:) = 0._wp      ! east
231         IF( nbondj == -1 .OR. nbondj == 2 )   ze3divh( :,2:nbghostcells+1        ,:) = 0._wp      ! south
232         IF( nbondj ==  1 .OR. nbondj == 2 )   ze3divh( :,nlcj-nbghostcells:nlcj-1,:) = 0._wp      ! north
233      ENDIF
234#endif
235      !
236      CALL lbc_lnk( ze3divh, 'T', 1. )
237      !
238      IF( ln_linssh ) THEN   ;   ik = 1   ! none zero velocity through the sea surface
239      ELSE                   ;   ik = 2   ! w=0 at the surface (set one for all in sbc_wave_init)
240      ENDIF
241      DO jk = jpkm1, ik, -1          ! integrate from the bottom the hor. divergence (NB: at k=jpk w is always zero)
242         wsd(:,:,jk) = wsd(:,:,jk+1) - ze3divh(:,:,jk)
243      END DO
244      !
245      IF( ln_bdy ) THEN
246         DO jk = 1, jpkm1
247            wsd(:,:,jk) = wsd(:,:,jk) * bdytmask(:,:)
248         END DO
249      ENDIF
250      !                       !==  Horizontal divergence of barotropic Stokes transport  ==!
251      div_sd(:,:) = 0._wp
252      DO jk = 1, jpkm1                                 !
253        div_sd(:,:) = div_sd(:,:) + ze3divh(:,:,jk)
254      END DO
255      !
256      CALL iom_put( "ustokes",  usd  )
257      CALL iom_put( "vstokes",  vsd  )
258      CALL iom_put( "wstokes",  wsd  )
259      !
260      DEALLOCATE( ze3divh )
261      DEALLOCATE( zk_t, zk_u, zk_v, zu0_sd, zv0_sd )
262      !
263   END SUBROUTINE sbc_stokes
264
265
266   SUBROUTINE sbc_wstress( )
267      !!---------------------------------------------------------------------
268      !!                     ***  ROUTINE sbc_wstress  ***
269      !!
270      !! ** Purpose :   Updates the ocean momentum modified by waves
271      !!
272      !! ** Method  : - Calculate u,v components of stress depending on stress
273      !!                model
274      !!              - Calculate the stress module
275      !!              - The wind module is not modified by waves
276      !! ** action 
277      !!---------------------------------------------------------------------
278      INTEGER  ::   jj, ji   ! dummy loop argument
279      !
280      IF( ln_tauwoc ) THEN
281         utau(:,:) = utau(:,:)*tauoc_wave(:,:)
282         vtau(:,:) = vtau(:,:)*tauoc_wave(:,:)
283         taum(:,:) = taum(:,:)*tauoc_wave(:,:)
284      ENDIF
285      !
286      IF( ln_tauw ) THEN
287         DO jj = 1, jpjm1
288            DO ji = 1, jpim1
289               ! Stress components at u- & v-points
290               utau(ji,jj) = 0.5_wp * ( tauw_x(ji,jj) + tauw_x(ji+1,jj) )
291               vtau(ji,jj) = 0.5_wp * ( tauw_y(ji,jj) + tauw_y(ji,jj+1) )
292               !
293               ! Stress module at t points
294               taum(ji,jj) = SQRT( tauw_x(ji,jj)*tauw_x(ji,jj) + tauw_y(ji,jj)*tauw_y(ji,jj) )
295            END DO
296         END DO
297         CALL lbc_lnk_multi( utau(:,:), 'U', -1. , vtau(:,:), 'V', -1. , taum(:,:) , 'T', -1. )
298      ENDIF
299      !
300   END SUBROUTINE sbc_wstress
301
302
303   SUBROUTINE sbc_wave( kt )
304      !!---------------------------------------------------------------------
305      !!                     ***  ROUTINE sbc_wave  ***
306      !!
307      !! ** Purpose :   read wave parameters from wave model  in netcdf files.
308      !!
309      !! ** Method  : - Read namelist namsbc_wave
310      !!              - Read Cd_n10 fields in netcdf files
311      !!              - Read stokes drift 2d in netcdf files
312      !!              - Read wave number in netcdf files
313      !!              - Compute 3d stokes drift using Breivik et al.,2014
314      !!                formulation
315      !! ** action 
316      !!---------------------------------------------------------------------
317      INTEGER, INTENT(in   ) ::   kt   ! ocean time step
318      !!---------------------------------------------------------------------
319      !
320      IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN     !==  Neutral drag coefficient  ==!
321         CALL fld_read( kt, nn_fsbc, sf_cd )             ! read from external forcing
322         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
323      ENDIF
324
325      IF( ln_tauwoc .AND. .NOT. cpl_tauwoc ) THEN  !==  Wave induced stress  ==!
326         CALL fld_read( kt, nn_fsbc, sf_tauwoc )         ! read wave norm stress from external forcing
327         tauoc_wave(:,:) = sf_tauwoc(1)%fnow(:,:,1)
328      ENDIF
329
330      IF( ln_tauw .AND. .NOT. cpl_tauw ) THEN      !==  Wave induced stress  ==!
331         CALL fld_read( kt, nn_fsbc, sf_tauw )           ! read ocean stress components from external forcing (T grid)
332         tauw_x(:,:) = sf_tauw(1)%fnow(:,:,1)
333         tauw_y(:,:) = sf_tauw(2)%fnow(:,:,1)
334      ENDIF
335
336      IF( ln_sdw )  THEN                           !==  Computation of the 3d Stokes Drift  ==!
337         !
338         IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled
339            CALL fld_read( kt, nn_fsbc, sf_sd )          ! read wave parameters from external forcing
340            IF( jp_hsw > 0 )   hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height
341            IF( jp_wmp > 0 )   wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period
342            IF( jp_wfr > 0 )   wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1)   ! Peak wave frequency
343            IF( jp_usd > 0 )   ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point
344            IF( jp_vsd > 0 )   vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point
345         ENDIF
346         !
347         ! Read also wave number if needed, so that it is available in coupling routines
348         IF( ln_zdfswm .AND. .NOT.cpl_wnum ) THEN
349            CALL fld_read( kt, nn_fsbc, sf_wn )          ! read wave parameters from external forcing
350            wnum(:,:) = sf_wn(1)%fnow(:,:,1)
351         ENDIF
352           
353         ! Calculate only if required fields have been read
354         ! In coupled wave model-NEMO case the call is done after coupling
355         !
356         IF( ( ll_st_bv_li   .AND. jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0 ) .OR. &
357           & ( ll_st_peakfr  .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0                ) ) CALL sbc_stokes()
358         !
359      ENDIF
360      !
361   END SUBROUTINE sbc_wave
362
363
364   SUBROUTINE sbc_wave_init
365      !!---------------------------------------------------------------------
366      !!                     ***  ROUTINE sbc_wave_init  ***
367      !!
368      !! ** Purpose :   read wave parameters from wave model  in netcdf files.
369      !!
370      !! ** Method  : - Read namelist namsbc_wave
371      !!              - Read Cd_n10 fields in netcdf files
372      !!              - Read stokes drift 2d in netcdf files
373      !!              - Read wave number in netcdf files
374      !!              - Compute 3d stokes drift using Breivik et al.,2014
375      !!                formulation
376      !! ** action 
377      !!---------------------------------------------------------------------
378      INTEGER ::   ierror, ios   ! local integer
379      INTEGER ::   ifpr
380      !!
381      CHARACTER(len=100)     ::  cn_dir                            ! Root directory for location of drag coefficient files
382      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i, slf_j     ! array of namelist informations on the fields to read
383      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  &
384                             &   sn_hsw, sn_wmp, sn_wfr, sn_wnum, &
385                             &   sn_tauwoc, sn_tauwx, sn_tauwy     ! informations about the fields to be read
386      !
387      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, &
388                             sn_wnum, sn_tauwoc, sn_tauwx, sn_tauwy
389      !!---------------------------------------------------------------------
390      !
391      REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model
392      READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901)
393901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp )
394         
395      REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model
396      READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 )
397902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp )
398      IF(lwm) WRITE ( numond, namsbc_wave )
399      !
400      IF( ln_cdgw ) THEN
401         IF( .NOT. cpl_wdrag ) THEN
402            ALLOCATE( sf_cd(1), STAT=ierror )               !* allocate and fill sf_wave with sn_cdg
403            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
404            !
405                                   ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   )
406            IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) )
407            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' )
408         ENDIF
409         ALLOCATE( cdn_wave(jpi,jpj) )
410      ENDIF
411
412      IF( ln_tauwoc ) THEN
413         IF( .NOT. cpl_tauwoc ) THEN
414            ALLOCATE( sf_tauwoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauwoc
415            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
416            !
417                                     ALLOCATE( sf_tauwoc(1)%fnow(jpi,jpj,1)   )
418            IF( sn_tauwoc%ln_tint )  ALLOCATE( sf_tauwoc(1)%fdta(jpi,jpj,1,2) )
419            CALL fld_fill( sf_tauwoc, (/ sn_tauwoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' )
420         ENDIF
421         ALLOCATE( tauoc_wave(jpi,jpj) )
422      ENDIF
423
424      IF( ln_tauw ) THEN
425         IF( .NOT. cpl_tauw ) THEN
426            ALLOCATE( sf_tauw(2), STAT=ierror )           !* allocate and fill sf_wave with sn_tauwx/y
427            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauw structure' )
428            !
429            ALLOCATE( slf_j(2) )
430            slf_j(1) = sn_tauwx
431            slf_j(2) = sn_tauwy
432                                    ALLOCATE( sf_tauw(1)%fnow(jpi,jpj,1)   )
433                                    ALLOCATE( sf_tauw(2)%fnow(jpi,jpj,1)   )
434            IF( slf_j(1)%ln_tint )  ALLOCATE( sf_tauw(1)%fdta(jpi,jpj,1,2) )
435            IF( slf_j(2)%ln_tint )  ALLOCATE( sf_tauw(2)%fdta(jpi,jpj,1,2) )
436            CALL fld_fill( sf_tauw, (/ slf_j /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
437         ENDIF
438         ALLOCATE( tauw_x(jpi,jpj) )
439         ALLOCATE( tauw_y(jpi,jpj) )
440      ENDIF
441
442      IF( ln_sdw ) THEN   ! Find out how many fields have to be read from file if not coupled
443         jpfld=0
444         jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0   ;   jp_wfr=0
445         IF( .NOT. cpl_sdrftx ) THEN
446            jpfld  = jpfld + 1
447            jp_usd = jpfld
448         ENDIF
449         IF( .NOT. cpl_sdrfty ) THEN
450            jpfld  = jpfld + 1
451            jp_vsd = jpfld
452         ENDIF
453         IF( .NOT. cpl_hsig  .AND. ll_st_bv_li  ) THEN
454            jpfld  = jpfld + 1
455            jp_hsw = jpfld
456         ENDIF
457         IF( .NOT. cpl_wper  .AND. ll_st_bv_li  ) THEN
458            jpfld  = jpfld + 1
459            jp_wmp = jpfld
460         ENDIF
461         IF( .NOT. cpl_wfreq .AND. ll_st_peakfr ) THEN
462            jpfld  = jpfld + 1
463            jp_wfr = jpfld
464         ENDIF
465
466         ! Read from file only the non-coupled fields
467         IF( jpfld > 0 ) THEN
468            ALLOCATE( slf_i(jpfld) )
469            IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd
470            IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd
471            IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw
472            IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp
473            IF( jp_wfr > 0 )   slf_i(jp_wfr) = sn_wfr
474
475            ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift
476            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
477            !
478            DO ifpr= 1, jpfld
479               ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) )
480               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) )
481            END DO
482            !
483            CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' )
484         ENDIF
485         ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) )
486         ALLOCATE( hsw  (jpi,jpj)    , wmp  (jpi,jpj)     )
487         ALLOCATE( wfreq(jpi,jpj) )
488         ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     )
489         ALLOCATE( div_sd(jpi,jpj) )
490         ALLOCATE( tsd2d (jpi,jpj) )
491
492         ut0sd(:,:) = 0._wp
493         vt0sd(:,:) = 0._wp
494         hsw(:,:) = 0._wp
495         wmp(:,:) = 0._wp
496
497         usd(:,:,:) = 0._wp
498         vsd(:,:,:) = 0._wp
499         wsd(:,:,:) = 0._wp
500         ! Wave number needed only if ln_zdfswm=T
501         IF( .NOT. cpl_wnum ) THEN
502            ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
503            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wave structure' )
504                                   ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   )
505            IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) )
506            CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
507         ENDIF
508         ALLOCATE( wnum(jpi,jpj) )
509      ENDIF
510      !
511   END SUBROUTINE sbc_wave_init
512
513   !!======================================================================
514END MODULE sbcwave
Note: See TracBrowser for help on using the repository browser.