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/2012/dev_r3389_INGV4_stokes/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2012/dev_r3389_INGV4_stokes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90 @ 3529

Last change on this file since 3529 was 3529, checked in by adani, 11 years ago

Read 2d Stokes Drift (U,V) and wavenumber, compute the 3D fields and the vertical component of the Stokes Drift.
File size: 10.4 KB
Line 
1MODULE sbcwave
2   !!======================================================================
3   !!                       ***  MODULE  sbcwave  ***
4   !! Wave module
5   !!======================================================================
6   !! History :  3.3.1  !   2011-09  (Adani M)  Original code: Drag Coefficient
7   !!         :  3.4    !   2012-10  (Adani M)                 Stokes Drift
8   !!----------------------------------------------------------------------
9   USE iom             ! I/O manager library
10   USE in_out_manager  ! I/O manager
11   USE lib_mpp         ! distribued memory computing library
12   USE fldread        ! read input fields
13   USE sbc_oce        ! Surface boundary condition: ocean fields
14
15   
16   !!----------------------------------------------------------------------
17   !!   sbc_wave       : read drag coefficient from wave model in netcdf files
18   !!----------------------------------------------------------------------
19
20   IMPLICIT NONE
21   PRIVATE
22
23   PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs
24   
25   INTEGER , PARAMETER ::   jpfld  = 3           ! maximum number of files to read for srokes drift
26   INTEGER , PARAMETER ::   jp_usd = 1           ! index of stokes drift  (i-component) (m/s)    at T-point
27   INTEGER , PARAMETER ::   jp_vsd = 2           ! index of stokes drift  (j-component) (m/s)    at T-point
28   INTEGER , PARAMETER ::   jp_wn  = 3           ! index of wave number                 (1/m)    at T-point
29   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient
30   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift
31   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave 
32   REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum 
33   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d 
34
35   !! * Substitutions
36#  include "domzgr_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
39   !! $Id: $
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE sbc_wave( kt )
45      !!---------------------------------------------------------------------
46      !!                     ***  ROUTINE sbc_apr  ***
47      !!
48      !! ** Purpose :   read drag coefficient from wave model  in netcdf files.
49      !!
50      !! ** Method  : - Read namelist namsbc_wave
51      !!              - Read Cd_n10 fields in netcdf files
52      !!              - Read stokes drift 2d in netcdf files
53      !!              - Read wave number      in netcdf files
54      !!              - Compute 3d stokes drift using monochromatic
55      !! ** action  :   
56      !!               
57      !!---------------------------------------------------------------------
58      USE oce,  ONLY : un,vn,hdivn,rotn
59      USE divcur
60      USE wrk_nemo
61#if defined key_bdy
62      USE bdy_oce, ONLY : bdytmask
63#endif
64      INTEGER, INTENT( in  ) ::  kt       ! ocean time step
65      INTEGER                ::  ierror   ! return error code
66      INTEGER                ::  ifpr, jj,ji,jk 
67      REAL(wp),DIMENSION(:,:,:),POINTER             ::  udummy,vdummy,hdivdummy,rotdummy
68      REAL                                          ::  z2dt,z1_2dt
69      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read
70      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files
71      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_wn   ! informations about the fields to be read
72      !!---------------------------------------------------------------------
73      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn
74      !!---------------------------------------------------------------------
75
76      !!----------------------------------------------------------------------
77      !
78      !
79      !                                         ! -------------------- !
80      IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
81         !                                      ! -------------------- !
82         !                                            !* set file information (default values)
83         ! ... default values (NB: frequency positive => hours, negative => months)
84         !              !   file   ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation !
85         !              !   name   !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    !
86         sn_cdg = FLD_N('cdg_wave'  ,    1     ,'drag_coeff',  .true.    , .false. ,   'daily'   , ''       , ''       )
87         sn_usd = FLD_N('sdw_wave'  ,    1     ,'u_sd2d',      .true.    , .false. ,   'daily'   , ''       , ''       )
88         sn_vsd = FLD_N('sdw_wave'  ,    1     ,'v_sd2d',      .true.    , .false. ,   'daily'   , ''       , ''       )
89         sn_wn = FLD_N( 'sdw_wave'  ,    1     ,'wave_num',    .true.    , .false. ,   'daily'   , ''       , ''       )
90         cn_dir = './'          ! directory in which the wave data are
91         
92
93         REWIND( numnam )                             !* read in namlist namsbc_wave
94         READ  ( numnam, namsbc_wave ) 
95         !
96
97         IF ( ln_cdgw ) THEN
98            ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
99            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
100            !
101                                   ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   )
102            IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) )
103            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
104            ALLOCATE( cdn_wave(jpi,jpj) )
105            cdn_wave(:,:) = 0.0
106        ENDIF
107         IF ( ln_sdw ) THEN
108            slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn
109            ALLOCATE( sf_sd(3), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
110            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
111            !
112            DO ifpr= 1, jpfld
113               ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) )
114               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) )
115            END DO
116            CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
117            ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj),uwavenum(jpi,jpj),vwavenum(jpi,jpj) )
118            ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) )
119            usd2d(:,:) = 0.0 ;  vsd2d(:,:) = 0.0 ; uwavenum(:,:) = 0.0 ; vwavenum(:,:) = 0.0
120            usd3d(:,:,:) = 0.0 ;vsd3d(:,:,:) = 0.0 ; wsd3d(:,:,:) = 0.0
121         ENDIF
122      ENDIF
123         !
124         !
125      IF ( ln_cdgw ) THEN
126         CALL fld_read( kt, nn_fsbc, sf_cd )      !* read drag coefficient from external forcing
127         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
128      ENDIF
129      IF ( ln_sdw )  THEN
130          CALL fld_read( kt, nn_fsbc, sf_sd )      !* read drag coefficient from external forcing
131
132         ! Interpolate wavenumber, stokes drift into the grid_V and grid_V
133         !-------------------------------------------------
134
135         DO jj = 1, jpjm1
136            DO ji = 1, jpim1
137               uwavenum(ji,jj)=0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) &
138               &                                + sf_sd(3)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) )
139
140               vwavenum(ji,jj)=0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) &
141               &                                + sf_sd(3)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) )
142
143               usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(1)%fnow(ji,jj,1) * tmask(ji,jj,1) &
144               &                                + sf_sd(1)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) )
145
146               vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(2)%fnow(ji,jj,1) * tmask(ji,jj,1) &
147               &                                + sf_sd(2)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) )
148            END DO
149         END DO
150
151          !Computation of the 3d Stokes Drift
152          DO jk = 1, jpk
153             DO jj = 1, jpj-1
154                DO ji = 1, jpi-1
155                   usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept(ji,jj,jk) , gdept(ji+1,jj  ,jk))))
156                   vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept(ji,jj,jk) , gdept(ji  ,jj+1,jk))))
157                END DO
158             END DO
159             usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept(jpi,:,jk)) )
160             vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept(:,jpj,jk)) )
161          END DO
162
163          CALL wrk_alloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
164         
165          udummy(:,:,:)=un(:,:,:)
166          vdummy(:,:,:)=vn(:,:,:)
167          hdivdummy(:,:,:)=hdivn(:,:,:)
168          rotdummy(:,:,:)=rotn(:,:,:)
169          un(:,:,:)=usd3d(:,:,:)
170          vn(:,:,:)=vsd3d(:,:,:)
171          CALL div_cur(kt)
172      !                                           !------------------------------!
173      !                                           !     Now Vertical Velocity    !
174      !                                           !------------------------------!
175          z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
176
177          z1_2dt = 1.e0 / z2dt
178          DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
179             ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
180             wsd3d(:,:,jk) = wsd3d(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
181                &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
182                &                         * tmask(:,:,jk) * z1_2dt
183#if defined key_bdy
184             wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:)
185#endif
186          END DO
187          hdivn(:,:,:)=hdivdummy(:,:,:)
188          rotn(:,:,:)=rotdummy(:,:,:)
189          vn(:,:,:)=vdummy(:,:,:)
190          un(:,:,:)=udummy(:,:,:)
191          CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy)
192      ENDIF
193   END SUBROUTINE sbc_wave
194     
195   !!======================================================================
196END MODULE sbcwave
Note: See TracBrowser for help on using the repository browser.