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.
dtauvd.F90 in branches/2013/dev_r3987_UKMO6_C1D/NEMOGCM/NEMO/OPA_SRC/C1D – NEMO

source: branches/2013/dev_r3987_UKMO6_C1D/NEMOGCM/NEMO/OPA_SRC/C1D/dtauvd.F90 @ 4144

Last change on this file since 4144 was 4144, checked in by rfurner, 10 years ago

Commit for 2013 changes; see #1085
File size: 13.3 KB
Line 
1MODULE dtauvd
2   !!======================================================================
3   !!                     ***  MODULE  dtauvd  ***
4   !! Ocean data  :  read ocean U & V current data from gridded data
5   !!======================================================================
6   !! History :  3.5   ! 2013-08  (D. Calvert)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   dta_uvd_init   : read namelist and allocate data structures
11   !!   dta_uvd        : read and time-interpolate ocean U & V current data
12   !!----------------------------------------------------------------------
13   USE oce             ! ocean dynamics and tracers
14   USE dom_oce         ! ocean space and time domain
15   USE fldread         ! read input fields
16   USE in_out_manager  ! I/O manager
17   USE phycst          ! physical constants
18   USE lib_mpp         ! MPP library
19   USE wrk_nemo        ! Memory allocation
20   USE timing          ! Timing
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   dta_uvd_init   ! called by nemogcm.F90
26   PUBLIC   dta_uvd        ! called by istate.F90 and dyndmp.90
27
28   LOGICAL , PUBLIC ::   ln_uvd_init   = .FALSE.      ! Flag to initialise with U & V current data
29   LOGICAL , PUBLIC ::   ln_uvd_dyndmp = .FALSE.      ! Flag for Newtonian damping toward U & V current data
30
31   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_uvd   ! structure for input U & V current (file information and data)
32
33   !! * Substitutions
34#  include "domzgr_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
37   !! $Id: dtauvd.F90 2392 2010-11-15 21:20:05Z gm $
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE dta_uvd_init( ld_dyndmp )
43      !!----------------------------------------------------------------------
44      !!                   ***  ROUTINE dta_uvd_init  ***
45      !!                   
46      !! ** Purpose :   initialization of U & V current input data
47      !!
48      !! ** Method  : - read namc1d_uvd namelist
49      !!              - allocate U & V current data structure
50      !!              - fld_fill data structure with namelist information
51      !!----------------------------------------------------------------------
52      LOGICAL, INTENT(in), OPTIONAL ::   ld_dyndmp         ! force the initialization when dyndmp is used
53      !
54      INTEGER ::   ierr0, ierr1, ierr2, ierr3              ! temporary integers
55      !
56      CHARACTER(len=100)            ::   cn_dir            ! Root directory for location of files to be used
57      TYPE(FLD_N), DIMENSION(2)     ::   suv_i             ! Combined U & V namelist information
58      TYPE(FLD_N)                   ::   sn_ucur, sn_vcur  ! U & V data namelist information
59      !!
60      NAMELIST/namc1d_uvd/ ln_uvd_init, ln_uvd_dyndmp, cn_dir, sn_ucur, sn_vcur
61      !!----------------------------------------------------------------------
62      !
63      IF( nn_timing == 1 )  CALL timing_start('dta_uvd_init')
64      !
65      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0
66      !
67      cn_dir = './'                 ! default: directory in which the model is executed
68      !
69      !                             !==   read namelist namc1d_uvd: flags and field info   ==!
70
71      !                             ! default values (NB: frequency positive => hours, negative => months)
72      !              !   file   ! frequency ! variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation !
73      !              !   name   !  (hours)  !  name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    !
74      sn_ucur = FLD_N( 'ucurrent',   -1     , 'u_current',  .false.  , .true.  ,  'monthly'  , ''       , 'Ume'    )
75      sn_vcur = FLD_N( 'vcurrent',   -1     , 'v_current',  .false.  , .true.  ,  'monthly'  , ''       , 'Vme'    )
76
77      REWIND( numnam )              ! Read in namelist namc1d_uvd
78      READ  ( numnam, namc1d_uvd )
79
80      !                             ! force the initialization when dyndmp is used
81      IF( PRESENT( ld_dyndmp ) )   ln_uvd_dyndmp = .TRUE.
82     
83      IF(lwp) THEN                  ! control print
84         WRITE(numout,*)
85         WRITE(numout,*) 'dta_uvd_init : U & V current data '
86         WRITE(numout,*) '~~~~~~~~~~~~ '
87         WRITE(numout,*) '   Namelist namc1d_uvd : Set flags'
88         WRITE(numout,*) '      Initialization of ocean U & V current with input data   ln_uvd_init   = ', ln_uvd_init
89         WRITE(numout,*) '      Damping of ocean U & V current toward input data        ln_uvd_dyndmp = ', ln_uvd_dyndmp
90         WRITE(numout,*)
91         IF( .NOT. ln_uvd_init .AND. .NOT. ln_uvd_dyndmp ) THEN
92            WRITE(numout,*)
93            WRITE(numout,*) '   U & V current data not used'
94         ENDIF
95      ENDIF
96      !                             ! no initialization when restarting
97      IF( ln_rstart .AND. ln_uvd_init ) THEN
98         CALL ctl_warn( 'dta_uvd_init: ocean restart and U & V current data initialization, ',   &
99            &           'we keep the restart U & V current values and set ln_uvd_init to FALSE' )
100         ln_uvd_init = .FALSE.
101      ENDIF
102
103      !
104      IF(  ln_uvd_init .OR. ln_uvd_dyndmp  ) THEN
105         !                          !==   allocate the data arrays   ==!
106         ALLOCATE( sf_uvd(2), STAT=ierr0 )
107         IF( ierr0 > 0 ) THEN
108            CALL ctl_stop( 'dta_uvd_init: unable to allocate sf_uvd structure' )             ;   RETURN
109         ENDIF
110         !
111                                 ALLOCATE( sf_uvd(1)%fnow(jpi,jpj,jpk)   , STAT=ierr0 )
112         IF( sn_ucur%ln_tint )   ALLOCATE( sf_uvd(1)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )
113                                 ALLOCATE( sf_uvd(2)%fnow(jpi,jpj,jpk)   , STAT=ierr2 )
114         IF( sn_vcur%ln_tint )   ALLOCATE( sf_uvd(2)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
115         !
116         IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN
117            CALL ctl_stop( 'dta_uvd_init : unable to allocate U & V current data arrays' )   ;   RETURN
118         ENDIF
119         !                          !==   fill sf_uvd with sn_ucur, sn_vcur and control print   ==!
120         suv_i(1) = sn_ucur   ;   suv_i(2) = sn_vcur
121         CALL fld_fill( sf_uvd, suv_i, cn_dir, 'dta_uvd', 'U & V current data', 'namc1d_uvd' )
122         !
123      ENDIF
124      !
125      IF( nn_timing == 1 )  CALL timing_stop('dta_uvd_init')
126      !
127   END SUBROUTINE dta_uvd_init
128
129
130   SUBROUTINE dta_uvd( kt, puvd )
131      !!----------------------------------------------------------------------
132      !!                   ***  ROUTINE dta_uvd  ***
133      !!                   
134      !! ** Purpose :   provides U & V current data at time step kt
135      !!
136      !! ** Method  : - call fldread routine
137      !!              - ORCA_R2: make some hand made alterations to the data (EMPTY)
138      !!              - s- or mixed s-zps coordinate: vertical interpolation onto model mesh
139      !!              - zps coordinate: vertical interpolation onto last partial level
140      !!              - ln_uvd_dyndmp=False: deallocate the U & V current data structure,
141      !!                                     as the data is no longer used
142      !!
143      !! ** Action  :   puvd,  U & V current data interpolated onto model mesh at time-step kt
144      !!----------------------------------------------------------------------
145      INTEGER                           , INTENT(in   ) ::   kt     ! ocean time-step
146      REAL(wp), DIMENSION(jpi,jpj,jpk,2), INTENT(  out) ::   puvd   ! U & V current data
147      !
148      INTEGER ::   ji, jj, jk, jl, jkk               ! dummy loop indicies
149      INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1  ! local integers
150      REAL(wp)::   zl, zi                            ! local floats
151      REAL(wp), POINTER, DIMENSION(:) ::  zup, zvp   ! 1D workspace
152      !!----------------------------------------------------------------------
153      !
154      IF( nn_timing == 1 )  CALL timing_start('dta_uvd')
155      !
156      CALL fld_read( kt, 1, sf_uvd )      !==   read U & V current data at time step kt   ==!
157      !
158      !
159      !                                   !==   ORCA_R2 configuration and U & V current damping   ==!
160      IF( cp_cfg == "orca" .AND. jp_cfg == 2 .AND. ln_uvd_dyndmp ) THEN    ! some hand made alterations
161         !!! EMPTY- to be added for running in 3D context !!!
162      ENDIF
163      !
164      puvd(:,:,:,1) = sf_uvd(1)%fnow(:,:,:)                 ! NO mask
165      puvd(:,:,:,2) = sf_uvd(2)%fnow(:,:,:) 
166      !
167      IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==!
168         !
169         CALL wrk_alloc( jpk, zup, zvp )
170         !
171         IF( kt == nit000 .AND. lwp )THEN
172            WRITE(numout,*)
173            WRITE(numout,*) 'dta_uvd: interpolate U & V current data onto the s- or mixed s-z-coordinate mesh'
174         ENDIF
175         !
176         DO jj = 1, jpj                   ! vertical interpolation of U & V current:
177            DO ji = 1, jpi                ! determines the interpolated U & V current profiles at each (i,j) point
178               DO jk = 1, jpk
179                  zl = fsdept_0(ji,jj,jk)
180                  IF    ( zl < gdept_0(1  ) ) THEN          ! extrapolate above the first level of data
181                     zup(jk) =  puvd(ji,jj,1    ,1)
182                     zvp(jk) =  puvd(ji,jj,1    ,2)
183                  ELSEIF( zl > gdept_0(jpk) ) THEN          ! extrapolate below the last level of data
184                     zup(jk) =  puvd(ji,jj,jpkm1,1)
185                     zvp(jk) =  puvd(ji,jj,jpkm1,2)
186                  ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
187                     DO jkk = 1, jpkm1                      ! when  gdept(jkk) < zl < gdept(jkk+1)
188                        IF( (zl-gdept_0(jkk)) * (zl-gdept_0(jkk+1)) <= 0._wp ) THEN
189                           zi = ( zl - gdept_0(jkk) ) / (gdept_0(jkk+1)-gdept_0(jkk))
190                           zup(jk) = puvd(ji,jj,jkk,1) + ( puvd(ji,jj,jkk+1,1 ) - puvd(ji,jj,jkk,1) ) * zi 
191                           zvp(jk) = puvd(ji,jj,jkk,2) + ( puvd(ji,jj,jkk+1,2 ) - puvd(ji,jj,jkk,2) ) * zi
192                        ENDIF
193                     END DO
194                  ENDIF
195               END DO
196               DO jk = 1, jpkm1           ! apply mask
197                  puvd(ji,jj,jk,1) = zup(jk) * umask(ji,jj,jk)
198                  puvd(ji,jj,jk,2) = zvp(jk) * vmask(ji,jj,jk)
199               END DO
200               puvd(ji,jj,jpk,1) = 0._wp
201               puvd(ji,jj,jpk,2) = 0._wp
202            END DO
203         END DO
204         !
205         CALL wrk_dealloc( jpk, zup, zvp )
206         !
207      ELSE                                !==   z- or zps- coordinate   ==!
208         !                             
209         puvd(:,:,:,1) = puvd(:,:,:,1) * umask(:,:,:)       ! apply mask
210         puvd(:,:,:,2) = puvd(:,:,:,2) * vmask(:,:,:)
211         !
212         IF( ln_zps ) THEN                ! zps-coordinate (partial steps) interpolation at the last ocean level
213            DO jj = 1, jpj
214               DO ji = 1, jpi
215                  ik = mbkt(ji,jj) 
216                  IF( ik > 1 ) THEN
217                     zl = ( gdept_0(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_0(ik) - gdept_0(ik-1) )
218                     puvd(ji,jj,ik,1) = (1.-zl) * puvd(ji,jj,ik,1) + zl * puvd(ji,jj,ik-1,1)
219                     puvd(ji,jj,ik,2) = (1.-zl) * puvd(ji,jj,ik,2) + zl * puvd(ji,jj,ik-1,2)
220                  ENDIF
221               END DO
222            END DO
223         ENDIF
224         !
225      ENDIF
226      !
227      IF( lwp .AND. kt == nit000 ) THEN   ! control print
228         WRITE(numout,*) ' U current '
229         WRITE(numout,*)
230         WRITE(numout,*)'  level = 1'
231         CALL prihre( puvd(:,:,1    ,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
232         WRITE(numout,*)'  level = ', jpk/2
233         CALL prihre( puvd(:,:,jpk/2,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
234         WRITE(numout,*)'  level = ', jpkm1
235         CALL prihre( puvd(:,:,jpkm1,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
236         WRITE(numout,*)
237         WRITE(numout,*) ' V current '
238         WRITE(numout,*)
239         WRITE(numout,*)'  level = 1'
240         CALL prihre( puvd(:,:,1    ,2), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
241         WRITE(numout,*)'  level = ', jpk/2
242         CALL prihre( puvd(:,:,jpk/2,2), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
243         WRITE(numout,*)'  level = ', jpkm1
244         CALL prihre( puvd(:,:,jpkm1,2), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
245         WRITE(numout,*)
246      ENDIF
247      !
248      IF( .NOT. ln_uvd_dyndmp    ) THEN   !==   deallocate U & V current structure   ==!
249         !                                !==   (data used only for initialization)  ==!
250         IF(lwp) WRITE(numout,*) 'dta_uvd: deallocate U & V current arrays as they are only used to initialize the run'
251                                   DEALLOCATE( sf_uvd(1)%fnow )     ! U current arrays in the structure
252         IF( sf_uvd(1)%ln_tint )   DEALLOCATE( sf_uvd(1)%fdta )
253                                   DEALLOCATE( sf_uvd(2)%fnow )     ! V current arrays in the structure
254         IF( sf_uvd(2)%ln_tint )   DEALLOCATE( sf_uvd(2)%fdta )
255                                   DEALLOCATE( sf_uvd         )     ! the structure itself
256      ENDIF
257      !
258      IF( nn_timing == 1 )  CALL timing_stop('dta_uvd')
259      !
260   END SUBROUTINE dta_uvd
261
262   !!======================================================================
263END MODULE dtauvd
Note: See TracBrowser for help on using the repository browser.