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/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90 @ 8553

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

Cap the value of tauoc as recieved via coupling or via forcing
File size: 27.4 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   !!             -   !   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 oce            ! ocean variables
19   USE sbc_oce        ! Surface boundary condition: ocean fields
20   USE bdy_oce        ! open boundary condition variables
21   USE domvvl         ! domain: variable volume layers
22   !
23   USE iom            ! I/O manager library
24   USE in_out_manager ! I/O manager
25   USE lib_mpp        ! distribued memory computing library
26   USE fldread        ! read input fields
27   USE wrk_nemo       !
28   USE phycst         ! physical constants
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   sbc_stokes      ! routine called in sbccpl
34   PUBLIC   sbc_wave        ! routine called in sbcmod
35   PUBLIC   sbc_wave_init   ! routine called in sbcmod
36   
37   ! Variables checking if the wave parameters are coupled (if not, they are read from file)
38   LOGICAL, PUBLIC ::   cpl_hsig   = .FALSE.
39   LOGICAL, PUBLIC ::   cpl_phioc  = .FALSE.
40   LOGICAL, PUBLIC ::   cpl_sdrft  = .FALSE.
41   LOGICAL, PUBLIC ::   cpl_wper   = .FALSE.
42   LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE.
43   LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE.
44   LOGICAL, PUBLIC ::   cpl_tauoc  = .FALSE.
45   LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE.
46
47   INTEGER ::   nn_sdrift      ! type of parameterization to calculate vertical Stokes drift
48   INTEGER, PARAMETER ::   jp_breivik  = 0     ! Breivik 2015: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]
49   INTEGER, PARAMETER ::   jp_phillips = 1     ! Phillips:     v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))]
50   INTEGER, PARAMETER ::   jp_peakph   = 2     ! Phillips using the peak wave number read from wave model instead of the inverse depth scale
51
52   INTEGER ::   jpfld    ! number of files to read for stokes drift
53   INTEGER ::   jp_usd   ! index of stokes drift  (i-component) (m/s)    at T-point
54   INTEGER ::   jp_vsd   ! index of stokes drift  (j-component) (m/s)    at T-point
55   INTEGER ::   jp_hsw   ! index of significant wave hight      (m)      at T-point
56   INTEGER ::   jp_wmp   ! index of mean wave period            (s)      at T-point
57   INTEGER ::   jp_wfr   ! index of wave peak frequency         (s^-1)   at T-point
58
59   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient
60   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift
61   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao
62   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean
63   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_phioc ! structure of input fields (file informations, fields read) wave to ocean energy
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(:,:)   ::   rn_crban            !: Craig and Banner constant for surface breaking waves mixing
68   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !:
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#  include "vectopt_loop_substitute.h90"
75   !!----------------------------------------------------------------------
76   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
77   !! $Id$
78   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
79   !!----------------------------------------------------------------------
80CONTAINS
81
82   SUBROUTINE sbc_stokes( )
83      !!---------------------------------------------------------------------
84      !!                     ***  ROUTINE sbc_stokes  ***
85      !!
86      !! ** Purpose :   compute the 3d Stokes Drift according to Breivik et al.,
87      !!                2014 (DOI: 10.1175/JPO-D-14-0020.1)
88      !!
89      !! ** Method  : - Calculate Stokes transport speed
90      !!              - Calculate horizontal divergence
91      !!              - Integrate the horizontal divergenze from the bottom
92      !! ** action 
93      !!---------------------------------------------------------------------
94      INTEGER  ::   jj, ji, jk   ! dummy loop argument
95      INTEGER  ::   ik           ! local integer
96      REAL(wp) ::  ztransp, zfac, ztemp
97      REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v, zda_u, zda_v
98      REAL(wp), DIMENSION(:,:)  , POINTER ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd   ! 2D workspace
99      REAL(wp), DIMENSION(:,:,:), POINTER ::   ze3divh                            ! 3D workspace
100
101      !!---------------------------------------------------------------------
102      !
103
104      CALL wrk_alloc( jpi,jpj,jpk,   ze3divh )
105      CALL wrk_alloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 
106      !
107      ! select parameterization for the calculation of vertical Stokes drift
108      ! exp. wave number at t-point
109      IF( nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips ) THEN   ! (Eq. (19) in Breivick et al. (2014) )
110         zfac = 2.0_wp * rpi / 16.0_wp
111         DO jj = 1, jpj               
112            DO ji = 1, jpi
113               ! Stokes drift velocity estimated from Hs and Tmean
114               ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp )
115               ! Stokes surface speed
116               tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj))
117               ! Wavenumber scale
118               zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp )
119            END DO
120         END DO
121         DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points
122            DO ji = 1, jpim1
123               zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
124               zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
125               !
126               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
127               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
128            END DO
129         END DO
130      ELSE IF( nn_sdrift==jp_peakph ) THEN    ! peak wave number calculated from the peak frequency received by the wave model
131         DO jj = 1, jpjm1             
132            DO ji = 1, jpim1
133               zk_u(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji+1,jj)*wfreq(ji+1,jj) ) / grav
134               zk_v(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji,jj+1)*wfreq(ji,jj+1) ) / grav
135               !
136               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
137               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
138            END DO
139         END DO
140      ENDIF
141      !
142      !                       !==  horizontal Stokes Drift 3D velocity  ==!
143      IF( nn_sdrift==jp_breivik ) THEN
144         DO jk = 1, jpkm1
145            DO jj = 2, jpjm1
146               DO ji = 2, jpim1
147                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )
148                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )
149                  !                         
150                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
151                  zkh_v = zk_v(ji,jj) * zdep_v
152                  !                                ! Depth attenuation
153                  zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u )
154                  zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v )
155                  !
156                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
157                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
158               END DO
159            END DO
160         END DO
161      ELSE IF( nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph ) THEN
162         DO jk = 1, jpkm1
163            DO jj = 2, jpjm1
164               DO ji = 2, jpim1
165                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )
166                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )
167                  !                         
168                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
169                  zkh_v = zk_v(ji,jj) * zdep_v
170                  !                                ! Depth attenuation
171                  zda_u = EXP( -2.0_wp*zkh_u ) - SQRT(2.0_wp*rpi*zkh_u) * ERFC(SQRT(2.0_wp*zkh_u))
172                  zda_v = EXP( -2.0_wp*zkh_v ) - SQRT(2.0_wp*rpi*zkh_v) * ERFC(SQRT(2.0_wp*zkh_v))
173                  !
174                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
175                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
176               END DO
177            END DO
178         END DO
179      ENDIF
180
181      CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. )
182      !
183      !                       !==  vertical Stokes Drift 3D velocity  ==!
184      !
185      DO jk = 1, jpkm1               ! Horizontal e3*divergence
186         DO jj = 2, jpj
187            DO ji = fs_2, jpi
188               ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * usd(ji,  jj,jk)    &
189                  &                 - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * usd(ji-1,jj,jk)    &
190                  &                 + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vsd(ji,jj  ,jk)    &
191                  &                 - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vsd(ji,jj-1,jk)  ) * r1_e12t(ji,jj)
192            END DO
193         END DO
194      END DO
195      !
196      IF( .NOT. AGRIF_Root() ) THEN
197         IF( nbondi ==  1 .OR. nbondi == 2 )   ze3divh(nlci-1,   :  ,:) = 0._wp      ! east
198         IF( nbondi == -1 .OR. nbondi == 2 )   ze3divh(  2   ,   :  ,:) = 0._wp      ! west
199         IF( nbondj ==  1 .OR. nbondj == 2 )   ze3divh(  :   ,nlcj-1,:) = 0._wp      ! north
200         IF( nbondj == -1 .OR. nbondj == 2 )   ze3divh(  :   ,  2   ,:) = 0._wp      ! south
201      ENDIF
202      !
203      CALL lbc_lnk( ze3divh, 'T', 1. )
204      !
205      IF( .NOT. lk_vvl ) THEN   ;   ik = 1   ! none zero velocity through the sea surface
206      ELSE                      ;   ik = 2   ! w=0 at the surface (set one for all in sbc_wave_init)
207      ENDIF
208      DO jk = jpkm1, ik, -1          ! integrate from the bottom the hor. divergence (NB: at k=jpk w is always zero)
209         wsd(:,:,jk) = wsd(:,:,jk+1) - ze3divh(:,:,jk)
210      END DO
211#if defined key_bdy
212      IF( lk_bdy ) THEN
213         DO jk = 1, jpkm1
214            wsd(:,:,jk) = wsd(:,:,jk) * bdytmask(:,:)
215         END DO
216      ENDIF
217#endif
218      !                       !==  Horizontal divergence of barotropic Stokes transport  ==!
219      div_sd(:,:) = 0._wp
220      DO jk = 1, jpkm1                                 !
221        div_sd(:,:) = div_sd(:,:) + ze3divh(:,:,jk)
222      END DO
223      !
224      CALL iom_put( "ustokes",  usd  )
225      CALL iom_put( "vstokes",  vsd  )
226      CALL iom_put( "wstokes",  wsd  )
227      !
228      CALL wrk_dealloc( jpi,jpj,jpk,   ze3divh )
229      CALL wrk_dealloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd )
230      !
231   END SUBROUTINE sbc_stokes
232
233
234   SUBROUTINE sbc_wave( kt )
235      !!---------------------------------------------------------------------
236      !!                     ***  ROUTINE sbc_wave  ***
237      !!
238      !! ** Purpose :   read wave parameters from wave model  in netcdf files.
239      !!
240      !! ** Method  : - Read namelist namsbc_wave
241      !!              - Read Cd_n10 fields in netcdf files
242      !!              - Read stokes drift 2d in netcdf files
243      !!              - Read wave number in netcdf files
244      !!              - Compute 3d stokes drift using Breivik et al.,2014
245      !!                formulation
246      !! ** action 
247      !!---------------------------------------------------------------------
248      INTEGER, INTENT(in   ) ::   kt   ! ocean time step
249      !!---------------------------------------------------------------------
250      !
251      IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN     !==  Neutral drag coefficient  ==!
252         CALL fld_read( kt, nn_fsbc, sf_cd )             ! read from external forcing
253         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
254         ! check that the drag coefficient contains proper information even if
255         ! the masks do not match - the momentum stress is not masked!
256         WHERE( cdn_wave < 0.0 ) cdn_wave = 1.5e-3
257         WHERE( cdn_wave > 1.0 ) cdn_wave = 1.5e-3
258      ENDIF
259
260      IF( ln_tauoc .AND. .NOT. cpl_tauoc ) THEN    !==  Wave induced stress  ==!
261         CALL fld_read( kt, nn_fsbc, sf_tauoc )          ! read wave norm stress from external forcing
262         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1)
263         WHERE( tauoc_wave <   0.0 ) tauoc_wave = 1.0
264         WHERE( tauoc_wave > 100.0 ) tauoc_wave = 1.0
265      ENDIF
266
267      IF( ln_phioc .AND. .NOT. cpl_phioc ) THEN    !==  Wave to ocean energy  ==!
268         CALL fld_read( kt, nn_fsbc, sf_phioc )          ! read wave to ocean energy from external forcing
269         rn_crban(:,:) = 29.0 * sf_phioc(1)%fnow(:,:,1)     ! ! Alfa is phioc*sqrt(rau0/zrhoa)  : rau0=water density, zhroa= air density
270         WHERE( rn_crban > 1.e8   ) rn_crban = 0.0    !remove first mask mistmatch points, then cap values in case of low friction velocity
271         WHERE( rn_crban < 0.0    ) rn_crban = 0.0
272         WHERE( rn_crban > 1000.0 ) rn_crban = 1000.0
273      ENDIF
274
275      IF( ln_sdw .OR. ln_rough )  THEN             !==  Computation of the 3d Stokes Drift  ==!
276         !
277         IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled
278            CALL fld_read( kt, nn_fsbc, sf_sd )          ! read wave parameters from external forcing
279            IF( jp_hsw > 0 ) THEN
280               hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height
281               WHERE( hsw > 100.0 ) hsw = 0.0
282               WHERE( hsw <   0.0 ) hsw = 0.0
283            ENDIF
284            IF( jp_wmp > 0 ) THEN
285               wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period
286               WHERE( wmp > 100.0 ) wmp = 0.0
287               WHERE( wmp <   0.0 ) wmp = 0.0
288            ENDIF
289            IF( jp_wfr > 0 ) THEN
290               wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1)   ! Peak wave frequency
291               WHERE( wfreq <    0.0 ) wfreq = 0.0 
292               WHERE( wfreq >  100.0 ) wfreq = 0.0
293            ENDIF
294            IF( jp_usd > 0 ) THEN
295               ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point
296               WHERE( ut0sd < -100.0 ) ut0sd = 1.0
297               WHERE( ut0sd >  100.0 ) ut0sd = 1.0
298            ENDIF
299            IF( jp_vsd > 0 ) THEN
300               vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point
301               WHERE( vt0sd < -100.0 ) vt0sd = 1.0
302               WHERE( vt0sd >  100.0 ) vt0sd = 1.0
303            ENDIF
304         ENDIF
305      ENDIF
306      !
307      IF( ln_sdw ) THEN
308         ! Read also wave number if needed, so that it is available in coupling routines
309         IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN
310            CALL fld_read( kt, nn_fsbc, sf_wn )          ! read wave parameters from external forcing
311            wnum(:,:) = sf_wn(1)%fnow(:,:,1)
312         ENDIF
313           
314         !                                         !==  Computation of the 3d Stokes Drift  ==!
315         !
316         IF( ((nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) .AND. &
317                          jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0) .OR. &
318             (nn_sdrift==jp_peakph .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0) ) &
319            CALL sbc_stokes()            ! Calculate only if required fields are read
320         !                               ! In coupled wave model-NEMO case the call is done after coupling
321         !
322      ENDIF
323      !
324   END SUBROUTINE sbc_wave
325
326
327   SUBROUTINE sbc_wave_init
328      !!---------------------------------------------------------------------
329      !!                     ***  ROUTINE sbc_wave_init  ***
330      !!
331      !! ** Purpose :   read wave parameters from wave model  in netcdf files.
332      !!
333      !! ** Method  : - Read namelist namsbc_wave
334      !!              - Read Cd_n10 fields in netcdf files
335      !!              - Read stokes drift 2d in netcdf files
336      !!              - Read wave number in netcdf files
337      !!              - Compute 3d stokes drift using Breivik et al.,2014
338      !!                formulation
339      !! ** action 
340      !!---------------------------------------------------------------------
341      INTEGER ::   ierror, ios   ! local integer
342      INTEGER ::   ifpr
343      !!
344      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files
345      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i     ! array of namelist informations on the fields to read
346      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_phioc, &
347                             &   sn_hsw, sn_wmp, sn_wfr, sn_wnum , &
348                             &   sn_tauoc      ! informations about the fields to be read
349      !
350      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, sn_wnum, sn_tauoc, sn_phioc, &
351                             ln_cdgw, ln_sdw, ln_stcor, ln_phioc, ln_tauoc, ln_zdfqiao, ln_rough,                 &
352                             nn_sdrift, nn_wmix
353      !!---------------------------------------------------------------------
354      !
355      REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model
356      READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901)
357901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp )
358         
359      REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model
360      READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 )
361902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp )
362      IF(lwm) WRITE ( numond, namsbc_wave )
363      !
364      IF(lwp) THEN               !* Parameter print
365         WRITE(numout,*)
366         WRITE(numout,*) 'sbc_wave_init: wave physics'
367         WRITE(numout,*) '~~~~~~~~'
368         WRITE(numout,*) '   Namelist namsbc_wave : set wave physics parameters'
369         WRITE(numout,*) '      Stokes drift corr. to vert. velocity    ln_sdw      = ', ln_sdw 
370         WRITE(numout,*) '        vertical parametrization              nn_sdrift   = ', nn_sdrift 
371         WRITE(numout,*) '      Stokes coriolis term                    ln_stcor    = ', ln_stcor 
372         WRITE(numout,*) '      wave modified ocean stress              ln_tauoc    = ', ln_tauoc 
373         WRITE(numout,*) '      wave to ocean energy                    ln_phioc    = ', ln_phioc
374         WRITE(numout,*) '        vertical mixing parametrization       nn_wmix     = ', nn_wmix 
375         WRITE(numout,*) '      neutral drag coefficient                ln_cdgw     = ', ln_cdgw
376         WRITE(numout,*) '      wave roughness length modification      ln_rough    = ', ln_rough 
377         WRITE(numout,*) '      Qiao vertical mixing formulation        ln_zdfqiao  = ', ln_zdfqiao
378      ENDIF
379
380      IF ( ln_wave ) THEN
381         ! Activated wave physics but no wave physics components activated
382         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_phioc .OR. ln_rough .OR. ln_zdfqiao) )   THEN 
383             CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F, ln_phioc=F ', &
384                                                      'ln_rough=F, ln_zdfqiao=F' ) 
385         ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN 
386             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
387         ENDIF
388         IF ( ln_cdgw .AND. .NOT.(nn_drag==jp_ukmo .OR. nn_drag==jp_std .OR. nn_drag==jp_const .OR. nn_drag==jp_mcore) ) & 
389             CALL ctl_stop( 'The chosen nn_drag for momentum calculation must be 0, 1, 2, or 3')
390         IF ( ln_cdgw .AND. ln_blk_core .AND. nn_drag==0 ) &
391             CALL ctl_stop( 'The chosen nn_drag for momentum calculation in core forcing must be 1, 2, or 3')
392         IF ( ln_cdgw .AND. ln_flx .AND. nn_drag==3 ) &
393             CALL ctl_stop( 'The chosen nn_drag for momentum calculation in direct forcing must be 0, 1, or 2')
394         IF( ln_phioc .AND. .NOT.(nn_wmix==jp_craigbanner .OR. nn_wmix==jp_janssen) ) & 
395            CALL ctl_stop( 'The chosen nn_wmix for wave vertical mixing must be 0, or 1' )
396         IF( ln_sdw .AND. .NOT.(nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph) ) & 
397            CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' )
398         IF( ln_zdfqiao .AND. .NOT.ln_sdw ) & 
399            CALL ctl_stop( 'Qiao vertical mixing can not be used without Stokes drift (ln_sdw)' )
400      ELSE
401         IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_phioc .OR. ln_rough .OR. ln_zdfqiao ) & 
402            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
403            &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
404            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 & 
405            &                  'or Stokes-Coriolis term (ln_stcor=T)',                         &
406            &                  'or ocean stress modification due to waves (ln_tauoc=T) ',      & 
407            &                  'or wave to ocean energy modification (ln_phioc=T) ',           & 
408            &                  'or wave surface roughness (ln_rough=T) ',                      & 
409            &                  'or Qiao vertical mixing formulation (ln_zdfqiao=T) ' )
410      ENDIF 
411      !
412      IF( ln_cdgw ) THEN
413         IF( .NOT. cpl_wdrag ) THEN
414            ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
415            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_cd structure' )
416            !
417                                   ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   )
418            IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) )
419            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
420         ENDIF
421         ALLOCATE( cdn_wave(jpi,jpj) )
422      ENDIF
423
424      IF( ln_tauoc ) THEN
425         IF( .NOT. cpl_tauoc ) THEN
426            ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc
427            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauoc structure' )
428            !
429                                    ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   )
430            IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) )
431            CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
432         ENDIF
433         ALLOCATE( tauoc_wave(jpi,jpj) )
434      ENDIF
435
436      IF( ln_phioc ) THEN
437         IF( .NOT. cpl_phioc ) THEN
438            ALLOCATE( sf_phioc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_phioc
439            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_phioc structure' )
440            !
441                                    ALLOCATE( sf_phioc(1)%fnow(jpi,jpj,1)   )
442            IF( sn_phioc%ln_tint )  ALLOCATE( sf_phioc(1)%fdta(jpi,jpj,1,2) )
443            CALL fld_fill( sf_phioc, (/ sn_phioc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
444         ENDIF
445         ALLOCATE( rn_crban(jpi,jpj) )
446      ENDIF
447
448      ! Find out how many fields have to be read from file if not coupled
449      jpfld=0
450      jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0   ;   jp_wfr=0
451      IF( ln_sdw ) THEN
452         IF( .NOT. cpl_sdrft ) THEN
453            jpfld  = jpfld + 1
454            jp_usd = jpfld
455            jpfld  = jpfld + 1
456            jp_vsd = jpfld
457         ENDIF
458         IF( .NOT. cpl_hsig .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN
459            jpfld  = jpfld + 1
460            jp_hsw = jpfld
461         ENDIF
462         IF( .NOT. cpl_wper .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN
463            jpfld  = jpfld + 1
464            jp_wmp = jpfld
465         ENDIF
466         IF( .NOT. cpl_wfreq .AND. nn_sdrift==jp_peakph ) THEN
467            jpfld  = jpfld + 1
468            jp_wfr = jpfld
469         ENDIF
470      ENDIF
471
472      IF( ln_rough .AND. .NOT. cpl_hsig .AND. jp_hsw==0 ) THEN
473         jpfld  = jpfld + 1
474         jp_hsw = jpfld
475      ENDIF
476
477      ! Read from file only the non-coupled fields
478      IF( jpfld > 0 ) THEN
479         ALLOCATE( slf_i(jpfld) )
480         IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd
481         IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd
482         IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw
483         IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp
484         IF( jp_wfr > 0 )   slf_i(jp_wfr) = sn_wfr
485         ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift
486         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_sd structure' )
487         !
488         DO ifpr= 1, jpfld
489            ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) )
490            IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) )
491         END DO
492         !
493         CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
494      ENDIF
495
496      IF( ln_sdw ) THEN
497         ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) )
498         ALLOCATE( wmp  (jpi,jpj)  )
499         ALLOCATE( wfreq (jpi,jpj) )
500         ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     )
501         ALLOCATE( div_sd(jpi,jpj) )
502         ALLOCATE( tsd2d (jpi,jpj) )
503         usd(:,:,:) = 0._wp
504         vsd(:,:,:) = 0._wp
505         wsd(:,:,:) = 0._wp
506         ! Wave number needed only if ln_zdfqiao=T
507         IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN
508            ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
509            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wn structure' )
510                                   ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   )
511            IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) )
512            CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'read wave input', 'namsbc_wave' )
513         ENDIF
514         ALLOCATE( wnum(jpi,jpj) )
515      ENDIF
516
517      IF( ln_sdw .OR. ln_rough ) THEN
518         ALLOCATE( hsw  (jpi,jpj) )
519      ENDIF
520      !
521   END SUBROUTINE sbc_wave_init
522
523   !!======================================================================
524END MODULE sbcwave
Note: See TracBrowser for help on using the repository browser.