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.
dtadyn.F90 in branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OFF_SRC – NEMO

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90 @ 2833

Last change on this file since 2833 was 2764, checked in by cetlod, 13 years ago

update the dtadyn.F90 routine, see ticket #823
  • Property svn:keywords set to Id
File size: 34.7 KB
Line 
1MODULE dtadyn
2   !!======================================================================
3   !!                       ***  MODULE  dtadyn  ***
4   !! Off-line : interpolation of the physical fields
5   !!======================================================================
6   !! History :   OPA  ! 1992-01 (M. Imbard) Original code
7   !!             8.0  ! 1998-04 (L.Bopp MA Foujols) slopes for isopyc.
8   !!              -   ! 1998-05 (L. Bopp) read output of coupled run
9   !!             8.2  ! 2001-01 (M. Levy et M. Benjelloul) add netcdf FORMAT
10   !!   NEMO      1.0  ! 2005-03 (O. Aumont and A. El Moussaoui) F90
11   !!              -   ! 2005-12 (C. Ethe) Adapted for DEGINT
12   !!             3.0  ! 2007-06 (C. Ethe) use of iom module
13   !!              -   ! 2007-09  (C. Ethe)  add swap_dyn_data
14   !!             3.3  ! 2010-11 (C. Ethe) Full reorganization of the off-line: phasing with the on-line
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   dta_dyn_init : initialization, namelist read, and parameters control
19   !!   dta_dyn      : Interpolation of the fields
20   !!----------------------------------------------------------------------
21   USE oce             ! ocean dynamics and tracers variables
22   USE c1d             ! 1D configuration: lk_c1d
23   USE dom_oce         ! ocean domain: variables
24   USE zdf_oce         ! ocean vertical physics: variables
25   USE sbc_oce         ! surface module: variables
26   USE phycst          ! physical constants
27   USE trabbl          ! active tracer: bottom boundary layer
28   USE ldfslp          ! lateral diffusion: iso-neutral slopes
29   USE ldfeiv          ! eddy induced velocity coef.
30   USE ldftra_oce      ! ocean tracer   lateral physics
31   USE zdfmxl          ! vertical physics: mixed layer depth
32   USE eosbn2          ! equation of state - Brunt Vaisala frequency
33   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
34   USE zpshde          ! z-coord. with partial steps: horizontal derivatives
35   USE in_out_manager  ! I/O manager
36   USE iom             ! I/O library
37   USE lib_mpp         ! distributed memory computing library
38   USE prtctl          !  print control
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   dta_dyn_init   ! called by opa.F90
44   PUBLIC   dta_dyn        ! called by step.F90
45
46   LOGICAL, PUBLIC ::   lperdyn = .TRUE.   !: boolean for periodic fields or not
47   LOGICAL, PUBLIC ::   lfirdyn = .TRUE.   !: boolean for the first call or not
48
49   INTEGER, PUBLIC ::   ndtadyn = 73       !: Number of dat in one year
50   INTEGER, PUBLIC ::   ndtatot = 73       !: Number of data in the input field
51   INTEGER, PUBLIC ::   nsptint = 1        !: type of spatial interpolation
52
53   CHARACTER(len=45) ::   cfile_grid_T = 'dyna_grid_T.nc'   ! name of the grid_T file
54   CHARACTER(len=45) ::   cfile_grid_U = 'dyna_grid_U.nc'   ! name of the grid_U file
55   CHARACTER(len=45) ::   cfile_grid_V = 'dyna_grid_V.nc'   ! name of the grid_V file
56   CHARACTER(len=45) ::   cfile_grid_W = 'dyna_grid_W.nc'   ! name of the grid_W file
57   
58   REAL(wp) ::   rnspdta    ! number of time step per 2 consecutives data
59   REAL(wp) ::   rnspdta2   ! rnspdta * 0.5
60
61   INTEGER ::   ndyn1, ndyn2    !
62   INTEGER ::   nlecoff = 0     ! switch for the first read
63   INTEGER ::   numfl_t, numfl_u, numfl_v, numfl_w
64
65   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: tdta       ! temperature at two consecutive times
66   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sdta       ! salinity at two consecutive times
67   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: udta       ! zonal velocity at two consecutive times
68   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vdta       ! meridional velocity at two consecutive times
69   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wdta       ! vertical velocity at two consecutive times
70   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: avtdta     ! vertical diffusivity coefficient
71
72   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: hmlddta    ! mixed layer depth at two consecutive times
73   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: wspddta    ! wind speed at two consecutive times
74   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: frlddta    ! sea-ice fraction at two consecutive times
75   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: empdta     ! E-P at two consecutive times
76   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: qsrdta     ! short wave heat flux at two consecutive times
77   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: bblxdta    ! bbl diffusive coef. in the x direction at 2 consecutive times
78   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: bblydta    ! bbl diffusive coef. in the y direction at 2 consecutive times
79   LOGICAL :: l_offbbl
80#if defined key_ldfslp && ! defined key_c1d
81   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta    ! zonal isopycnal slopes
82   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta    ! meridional isopycnal slopes
83   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta   ! zonal diapycnal slopes
84   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta   ! meridional diapycnal slopes
85#endif
86#if ! defined key_degrad &&  defined key_traldf_c2d && defined key_traldf_eiv
87   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: aeiwdta    ! G&M coefficient
88#endif
89#if defined key_degrad
90   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: ahtudta, ahtvdta, ahtwdta   ! Lateral diffusivity
91# if defined key_traldf_eiv
92   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: aeiudta, aeivdta, aeiwdta   ! G&M coefficient
93# endif
94#endif
95
96   !! * Substitutions
97#  include "domzgr_substitute.h90"
98#  include "vectopt_loop_substitute.h90"
99   !!----------------------------------------------------------------------
100   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
101   !! $Id$
102   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
103   !!----------------------------------------------------------------------
104CONTAINS
105
106   SUBROUTINE dta_dyn( kt )
107      !!----------------------------------------------------------------------
108      !!                  ***  ROUTINE dta_dyn  ***
109      !!
110      !! ** Purpose :   Prepares dynamics and physics fields from an NEMO run
111      !!              for an off-line simulation of passive tracers
112      !!
113      !! ** Method : calculates the position of DATA to read READ DATA
114      !!             (example month changement) computes slopes IF needed
115      !!             interpolates DATA IF needed
116      !!----------------------------------------------------------------------
117      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
118      !!
119      INTEGER  ::   iper, iperm1, iswap, izt   ! local integers
120      REAL(wp) ::   zt, zweigh                 ! local scalars
121      !!----------------------------------------------------------------------
122
123      zt     = ( REAL(kt,wp) + rnspdta2 ) / rnspdta
124      izt    = INT( zt )
125      zweigh = zt - REAL( INT(zt), wp )
126
127      IF( lperdyn ) THEN   ;   iperm1 = MOD( izt, ndtadyn )
128      ELSE                 ;   iperm1 = MOD( izt, ndtatot - 1 ) + 1
129      ENDIF
130
131      iper = iperm1 + 1
132      IF( iperm1 == 0 ) THEN
133          IF( lperdyn ) THEN
134              iperm1 = ndtadyn
135          ELSE
136              IF( lfirdyn ) THEN
137                  IF(lwp) WRITE (numout,*) 'dta_dyn:  dynamic file is not periodic with or without interpolation    &
138                     &                                we take the first value for the previous period iperm1 = 0  '
139              END IF
140          END IF
141      END IF
142
143      iswap  = 0
144
145      ! 1. First call lfirdyn = true
146      ! ----------------------------
147
148      IF( lfirdyn ) THEN
149         ndyn1 = iperm1         ! store the information of the period read
150         ndyn2 = iper
151         
152         IF(lwp) THEN
153            WRITE (numout,*) ' dynamics data read for the period ndyn1 =', ndyn1,   &
154               &             ' and for the period ndyn2 = ', ndyn2
155            WRITE (numout,*) ' time step is : ', kt
156            WRITE (numout,*) ' we have ndtadyn = ', ndtadyn, ' records in the dynamic file for one year'
157         END IF
158         !
159         CALL dynrea( kt, MAX( 1, iperm1) )           ! data read for the iperm1 period
160         !
161         CALL swap_dyn_data            ! swap from record 2 to 1
162         !
163         iswap = 1        !  indicates swap
164         !
165         CALL dynrea( kt, iper )       ! data read for the iper period
166         !
167         lfirdyn = .FALSE.    ! trace the first call
168      ENDIF
169      !
170      ! And now what we have to do at every time step
171      ! check the validity of the period in memory
172      !
173      IF( iperm1 /= ndyn1 ) THEN 
174         !
175         IF( iperm1 == 0 ) THEN
176            IF(lwp) THEN
177               WRITE (numout,*) ' dynamic file is not periodic with periodic interpolation'
178               WRITE (numout,*) ' we take the last value for the last period '
179               WRITE (numout,*) ' iperm1 = 12,   iper = 13  '
180            ENDIF
181            iperm1 = 12
182            iper   = 13
183         ENDIF
184         !
185         CALL swap_dyn_data         ! We have to prepare a new read of data : swap from record 2 to 1
186         !
187         iswap = 1                  !  indicates swap
188         !
189         CALL dynrea( kt, iper )    ! data read for the iper period
190         !
191         ndyn1 = ndyn2         ! store the information of the period read
192         ndyn2 = iper
193         !
194         IF(lwp) THEN
195            WRITE (numout,*) ' dynamics data read for the period ndyn1 =', ndyn1,   &
196               &             ' and for the period ndyn2 = ', ndyn2
197            WRITE (numout,*) ' time step is : ', kt
198         END IF
199         !
200      END IF
201      !
202      ! Compute the data at the given time step
203      !----------------------------------------     
204
205      IF( nsptint == 0 ) THEN          ! No space interpolation, data are probably correct
206         !                             ! We have to initialize data if we have changed the period         
207         CALL assign_dyn_data
208      ELSEIF( nsptint == 1 ) THEN      ! linear interpolation
209         CALL linear_interp_dyn_data( zweigh )
210      ELSE                             ! other interpolation
211         WRITE (numout,*) ' this kind of interpolation do not exist at the moment : we stop'
212         STOP 'dtadyn'         
213      END IF
214      !
215      CALL eos( tsn, rhd, rhop )       ! In any case, we need rhop
216      !
217#if ! defined key_degrad && defined key_traldf_c2d
218      !                                ! In case of 2D varying coefficients, we need aeiv and aeiu
219      IF( lk_traldf_eiv )   CALL dta_eiv( kt )      ! eddy induced velocity coefficient
220#endif
221      !
222      IF( .NOT. l_offbbl ) THEN       ! Compute bbl coefficients if needed
223         tsb(:,:,:,:) = tsn(:,:,:,:)
224         CALL bbl( kt, 'TRC')
225      END IF
226      !
227      IF(ln_ctl) THEN
228         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
229         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
230         CALL prt_ctl(tab3d_1=un               , clinfo1=' un      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
231         CALL prt_ctl(tab3d_1=vn               , clinfo1=' vn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
232         CALL prt_ctl(tab3d_1=wn               , clinfo1=' wn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
233         CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
234         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 )
235         CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask, ovlap=1 )
236         CALL prt_ctl(tab2d_1=emps             , clinfo1=' emps    - : ', mask1=tmask, ovlap=1 )
237         CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask, ovlap=1 )
238         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask, ovlap=1 )
239      ENDIF
240      !
241   END SUBROUTINE dta_dyn
242
243
244   INTEGER FUNCTION dta_dyn_alloc()
245      !!---------------------------------------------------------------------
246      !!                 ***  ROUTINE dta_dyn_alloc  ***
247      !!---------------------------------------------------------------------
248
249      ALLOCATE( tdta    (jpi,jpj,jpk,2), sdta    (jpi,jpj,jpk,2),    &
250         &      udta    (jpi,jpj,jpk,2), vdta    (jpi,jpj,jpk,2),    &
251         &      wdta    (jpi,jpj,jpk,2), avtdta  (jpi,jpj,jpk,2),    &
252#if defined key_ldfslp && ! defined key_c1d
253         &      uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    &
254         &      wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2),    &
255#endif
256#if defined key_degrad
257         &      ahtudta (jpi,jpj,jpk,2), ahtvdta (jpi,jpj,jpk,2),    &
258         &      ahtwdta (jpi,jpj,jpk,2),                             &
259# if defined key_traldf_eiv
260         &      aeiudta (jpi,jpj,jpk,2), aeivdta (jpi,jpj,jpk,2),    &
261         &      aeiwdta (jpi,jpj,jpk,2),                             &
262# endif
263#endif
264#if ! defined key_degrad &&  defined key_traldf_c2d && defined key_traldf_eiv
265         &      aeiwdta (jpi,jpj,    2),                             &
266#endif
267         &      hmlddta (jpi,jpj,    2), wspddta (jpi,jpj,    2),    &
268         &      frlddta (jpi,jpj,    2), qsrdta  (jpi,jpj,    2),    &
269         &      empdta  (jpi,jpj,    2),                         STAT=dta_dyn_alloc ) 
270         !
271      IF( dta_dyn_alloc /= 0 )   CALL ctl_warn('dta_dyn_alloc: failed to allocate facvol array')
272      !
273   END FUNCTION dta_dyn_alloc
274
275
276   SUBROUTINE dynrea( kt, kenr )
277      !!----------------------------------------------------------------------
278      !!                  ***  ROUTINE dynrea  ***
279      !!
280      !! ** Purpose : READ dynamics fiels from OPA9 netcdf output
281      !!
282      !! ** Method : READ the kenr records of DATA and store in udta(...,2), .... 
283      !!----------------------------------------------------------------------
284      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
285      USE wrk_nemo, ONLY: zu      => wrk_3d_3  , zv    => wrk_3d_4 , zw   => wrk_3d_5
286      USE wrk_nemo, ONLY: zt      => wrk_3d_6  , zs    => wrk_3d_7 , zavt => wrk_3d_8 
287      USE wrk_nemo, ONLY: zemp    => wrk_2d_11 , zqsr  => wrk_2d_12, zmld => wrk_2d_13
288      USE wrk_nemo, ONLY: zice    => wrk_2d_14 , zwspd => wrk_2d_15 
289      USE wrk_nemo, ONLY: ztaux   => wrk_2d_16 , ztauy => wrk_2d_17
290      USE wrk_nemo, ONLY: zbblx   => wrk_2d_18 , zbbly => wrk_2d_19
291      USE wrk_nemo, ONLY: zaeiw2d => wrk_2d_10
292      USE wrk_nemo, ONLY: ztsn    => wrk_4d_1
293      !
294      INTEGER, INTENT(in) ::   kt, kenr   ! time index
295      !!
296      INTEGER ::  jkenr
297#if defined key_degrad
298      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zahtu, zahtv, zahtw   ! Lateral diffusivity
299# if defined key_traldf_eiv
300      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zaeiu, zaeiv, zaeiw   ! G&M coefficient
301# endif
302#endif
303      !!----------------------------------------------------------------------
304      !
305      IF( wrk_in_use(3, 3,4,5,6,7,8) .OR. &
306          wrk_in_use(4, 1)                             .OR. &
307          wrk_in_use(2, 10,11,12,13,14,15,16,17,18,19)               ) THEN
308         CALL ctl_stop('domrea/dta_dyn: requested workspace arrays unavailable')   ;   RETURN
309      ENDIF
310
311#if defined key_degrad
312      ALLOCATE( zahtu(jpi,jpj,jpk), zahtv(jpi,jpj,jpk), zahtw(jpi,jpj,jpk) ) 
313# if defined key_traldf_eiv
314      ALLOCATE( zaeiu(jpi,jpj,jpk), zaeiv(jpi,jpj,jpk), zaeiw(jpi,jpj,jpk) )
315# endif
316#endif
317     
318      ! cas d'un fichier non periodique : on utilise deux fois le premier et
319      ! le dernier champ temporel
320
321      jkenr = kenr
322
323      IF(lwp) THEN
324         WRITE(numout,*)
325         WRITE(numout,*) 'Dynrea : read dynamical fields, kenr = ', jkenr
326         WRITE(numout,*) '~~~~~~~'
327#if defined key_degrad
328         WRITE(numout,*) ' Degraded fields'
329#endif
330         WRITE(numout,*)
331      ENDIF
332
333
334      IF( kt == nit000 .AND. nlecoff == 0 ) THEN
335         nlecoff = 1
336         CALL  iom_open ( cfile_grid_T, numfl_t )
337         CALL  iom_open ( cfile_grid_U, numfl_u )
338         CALL  iom_open ( cfile_grid_V, numfl_v )
339         CALL  iom_open ( cfile_grid_W, numfl_w )
340      ENDIF
341
342      ! file grid-T
343      !---------------
344      CALL iom_get( numfl_t, jpdom_data, 'votemper', zt   (:,:,:), jkenr )
345      CALL iom_get( numfl_t, jpdom_data, 'vosaline', zs   (:,:,:), jkenr )
346      CALL iom_get( numfl_t, jpdom_data, 'somixhgt', zmld (:,:  ), jkenr )
347      CALL iom_get( numfl_t, jpdom_data, 'sowaflcd', zemp (:,:  ), jkenr )
348      CALL iom_get( numfl_t, jpdom_data, 'soshfldo', zqsr (:,:  ), jkenr )
349      CALL iom_get( numfl_t, jpdom_data, 'soicecov', zice (:,:  ), jkenr )
350      IF( iom_varid( numfl_t, 'sowindsp', ldstop = .FALSE. ) > 0 ) THEN
351         CALL iom_get( numfl_t, jpdom_data, 'sowindsp', zwspd(:,:), jkenr ) 
352      ELSE
353         CALL iom_get( numfl_u, jpdom_data, 'sozotaux', ztaux(:,:), jkenr )
354         CALL iom_get( numfl_v, jpdom_data, 'sometauy', ztauy(:,:), jkenr )
355         CALL tau2wnd( ztaux, ztauy, zwspd )
356      ENDIF
357      ! files grid-U / grid_V
358      CALL iom_get( numfl_u, jpdom_data, 'vozocrtx', zu   (:,:,:), jkenr )
359      CALL iom_get( numfl_v, jpdom_data, 'vomecrty', zv   (:,:,:), jkenr )
360#if defined key_trabbl
361      IF( .NOT. lk_c1d .AND. nn_bbl_ldf == 1 ) THEN
362         IF( iom_varid( numfl_u, 'ahu_bbl', ldstop = .FALSE. ) > 0  .AND. &
363         &   iom_varid( numfl_v, 'ahv_bbl', ldstop = .FALSE. ) > 0 ) THEN
364             CALL iom_get( numfl_u, jpdom_data, 'ahu_bbl', zbblx(:,:), jkenr )
365             CALL iom_get( numfl_v, jpdom_data, 'ahv_bbl', zbbly(:,:), jkenr )
366             l_offbbl = .TRUE.
367         ENDIF
368      ENDIF
369#endif 
370
371      ! file grid-W
372      ! CALL iom_get ( numfl_w, jpdom_data, 'vovecrtz', zw   (:,:,:), jkenr )
373      ! Computation of vertical velocity using horizontal divergence
374      CALL wzv( zu, zv, zw )
375
376      IF( iom_varid( numfl_w, 'voddmavs', ldstop = .FALSE. ) > 0 ) THEN          ! avs exist: it is used
377         CALL iom_get( numfl_w, jpdom_data, 'voddmavs', zavt (:,:,:), jkenr )
378      ELSE                                                                       ! no avs: use avt
379         CALL iom_get( numfl_w, jpdom_data, 'votkeavt', zavt (:,:,:), jkenr )
380      ENDIF
381
382#if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv
383      CALL iom_get( numfl_w, jpdom_data, 'soleaeiw', zaeiw2d(:,: ), jkenr )
384#endif
385
386#if defined key_degrad
387      CALL iom_get( numfl_u, jpdom_data, 'vozoahtu', zahtu(:,:,:), jkenr )
388      CALL iom_get( numfl_v, jpdom_data, 'vomeahtv', zahtv(:,:,:), jkenr )
389      CALL iom_get( numfl_w, jpdom_data, 'voveahtw', zahtw(:,:,:), jkenr )
390#  if defined key_traldf_eiv
391      CALL iom_get( numfl_u, jpdom_data, 'vozoaeiu', zaeiu(:,:,:), jkenr )
392      CALL iom_get( numfl_v, jpdom_data, 'vomeaeiv', zaeiv(:,:,:), jkenr )
393      CALL iom_get( numfl_w, jpdom_data, 'voveaeiw', zaeiw(:,:,:), jkenr )
394#  endif
395#endif
396
397      udta  (:,:,:,2) = zu  (:,:,:) * umask(:,:,:)
398      vdta  (:,:,:,2) = zv  (:,:,:) * vmask(:,:,:) 
399      wdta  (:,:,:,2) = zw  (:,:,:) * tmask(:,:,:)
400      tdta  (:,:,:,2) = zt  (:,:,:) * tmask(:,:,:)
401      sdta  (:,:,:,2) = zs  (:,:,:) * tmask(:,:,:)
402      avtdta(:,:,:,2) = zavt(:,:,:) * tmask(:,:,:)
403
404#if defined key_ldfslp && ! defined key_c1d
405      ! Computes slopes (here tsn and avt are used as workspace)
406      ztsn (:,:,:,jp_tem) = tdta  (:,:,:,2)
407      ztsn (:,:,:,jp_sal) = sdta  (:,:,:,2)
408      avt(:,:,:)          = avtdta(:,:,:,2)
409     
410      CALL eos( ztsn, rhd, rhop )   ! Time-filtered in situ density
411      CALL bn2( ztsn, rn2 )         ! before Brunt-Vaisala frequency
412      IF( ln_zps )   &
413         &   CALL zps_hde( kt, jpts, ztsn, gtsu, gtsv,  &  ! Partial steps: before Horizontal DErivative
414         &                           rhd, gru , grv   )    ! of t, s, rd at the bottom ocean level
415      CALL zdf_mxl( kt )           ! mixed layer depth
416      CALL ldf_slp( kt, rhd, rn2 )
417         
418      uslpdta (:,:,:,2) = uslp (:,:,:)
419      vslpdta (:,:,:,2) = vslp (:,:,:)
420      wslpidta(:,:,:,2) = wslpi(:,:,:)
421      wslpjdta(:,:,:,2) = wslpj(:,:,:)
422#endif
423
424#if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv
425      aeiwdta(:,:,2)  = zaeiw2d(:,:) * tmask(:,:,1)
426#endif
427
428#if defined key_degrad
429        ahtudta(:,:,:,2) = zahtu(:,:,:) * umask(:,:,:)
430        ahtvdta(:,:,:,2) = zahtv(:,:,:) * vmask(:,:,:)
431        ahtwdta(:,:,:,2) = zahtw(:,:,:) * tmask(:,:,:)
432#  if defined key_traldf_eiv
433        aeiudta(:,:,:,2) = zaeiu(:,:,:) * umask(:,:,:)
434        aeivdta(:,:,:,2) = zaeiv(:,:,:) * vmask(:,:,:)
435        aeiwdta(:,:,:,2) = zaeiw(:,:,:) * tmask(:,:,:)
436#  endif
437#endif
438
439      ! fluxes
440      !
441      wspddta(:,:,2)  = zwspd(:,:) * tmask(:,:,1)
442      frlddta(:,:,2)  = zice (:,:) * tmask(:,:,1)
443      empdta (:,:,2)  = zemp (:,:) * tmask(:,:,1)
444      qsrdta (:,:,2)  = zqsr (:,:) * tmask(:,:,1)
445      hmlddta(:,:,2)  = zmld (:,:) * tmask(:,:,1)
446
447#if defined key_trabbl
448      IF( l_offbbl ) THEN
449         bblxdta(:,:,2) = zbblx(:,:)  * umask(:,:,1)
450         bblydta(:,:,2) = zbbly(:,:)  * vmask(:,:,1)
451      ENDIF
452#endif
453     
454      IF( kt == nitend ) THEN
455         CALL iom_close ( numfl_t )
456         CALL iom_close ( numfl_u )
457         CALL iom_close ( numfl_v )
458         CALL iom_close ( numfl_w )
459      ENDIF
460      !     
461      IF( wrk_not_released(3, 3,4,5,6,7,8) .OR. &
462          wrk_not_released(4, 1                            ) .OR. &
463          wrk_not_released(2, 10,11,12,13,14,15,16,17,18,19)                ) THEN
464         CALL ctl_stop('domrea/dta_dyn: failed to release workspace arrays')
465      END IF
466#if defined key_degrad
467      DEALLOCATE( zahtu )   ;   DEALLOCATE( zahtv )   ;   DEALLOCATE( zahtw )
468# if defined key_traldf_eiv
469      DEALLOCATE( zaeiu )   ;   DEALLOCATE( zaeiv )   ;   DEALLOCATE( zaeiw )
470# endif
471#endif
472      !
473   END SUBROUTINE dynrea
474
475
476   SUBROUTINE dta_dyn_init
477      !!----------------------------------------------------------------------
478      !!                  ***  ROUTINE dta_dyn_init  ***
479      !!
480      !! ** Purpose :   initializations of parameters for the interpolation
481      !!
482      !! ** Method :
483      !!----------------------------------------------------------------------
484      REAL(wp) :: znspyr   !: number of time step per year
485      !
486      NAMELIST/namdyn/ ndtadyn, ndtatot, nsptint, lperdyn,  &
487         &             cfile_grid_T, cfile_grid_U, cfile_grid_V, cfile_grid_W
488      !!----------------------------------------------------------------------
489      !
490      IF( dta_dyn_alloc() /= 0 )  CALL ctl_stop( 'STOP', 'dta_dyn_alloc: unable to allocate standard ocean arrays' )
491      !
492      REWIND( numnam )              ! Read Namelist namdyn : Lateral physics on tracers
493      READ  ( numnam, namdyn )
494      !
495      IF(lwp) THEN                  ! control print
496         WRITE(numout,*)
497         WRITE(numout,*) 'namdyn : offline dynamical selection'
498         WRITE(numout,*) '~~~~~~~'
499         WRITE(numout,*) '  Namelist namdyn : set parameters for the lecture of the dynamical fields'
500         WRITE(numout,*) 
501         WRITE(numout,*) ' number of elements in the FILE for a year  ndtadyn = ' , ndtadyn
502         WRITE(numout,*) ' total number of elements in the FILE       ndtatot = ' , ndtatot
503         WRITE(numout,*) ' type of interpolation                      nsptint = ' , nsptint
504         WRITE(numout,*) ' loop on the same FILE                      lperdyn = ' , lperdyn
505         WRITE(numout,*) '  '
506         WRITE(numout,*) ' name of grid_T file                   cfile_grid_T = ', TRIM(cfile_grid_T)   
507         WRITE(numout,*) ' name of grid_U file                   cfile_grid_U = ', TRIM(cfile_grid_U) 
508         WRITE(numout,*) ' name of grid_V file                   cfile_grid_V = ', TRIM(cfile_grid_V) 
509         WRITE(numout,*) ' name of grid_W file                   cfile_grid_W = ', TRIM(cfile_grid_W)     
510         WRITE(numout,*) ' '
511      ENDIF
512      !
513      znspyr   = nyear_len(1) * rday / rdt 
514      rnspdta  = znspyr / REAL( ndtadyn, wp )
515      rnspdta2 = rnspdta * 0.5 
516      !
517      CALL dta_dyn( nit000 )
518      !
519   END SUBROUTINE dta_dyn_init
520
521
522   SUBROUTINE wzv( pu, pv, pw )
523      !!----------------------------------------------------------------------
524      !!                    ***  ROUTINE wzv  ***
525      !!
526      !! ** Purpose :   Compute the now vertical velocity after the array swap
527      !!
528      !! ** Method  : - compute the now divergence given by :
529      !!         * z-coordinate ONLY !!!!
530      !!         hdiv = 1/(e1t*e2t) [ di(e2u  u) + dj(e1v  v) ]
531      !!     - Using the incompressibility hypothesis, the vertical
532      !!      velocity is computed by integrating the horizontal divergence
533      !!      from the bottom to the surface.
534      !!        The boundary conditions are w=0 at the bottom (no flux).
535      !!----------------------------------------------------------------------
536      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pu, pv    !:  horizontal velocities
537      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) :: pw        !:  verticla velocity
538      !!
539      INTEGER  ::  ji, jj, jk
540      REAL(wp) ::  zu, zu1, zv, zv1, zet
541      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhdiv     !:  horizontal divergence
542      !!----------------------------------------------------------------------
543      !
544      ! Computation of vertical velocity using horizontal divergence
545      zhdiv(:,:,:) = 0.
546      DO jk = 1, jpkm1
547         DO jj = 2, jpjm1
548            DO ji = fs_2, fs_jpim1   ! vector opt.
549               zu  = pu(ji  ,jj  ,jk) * umask(ji  ,jj  ,jk) * e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk)
550               zu1 = pu(ji-1,jj  ,jk) * umask(ji-1,jj  ,jk) * e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk)
551               zv  = pv(ji  ,jj  ,jk) * vmask(ji  ,jj  ,jk) * e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk)
552               zv1 = pv(ji  ,jj-1,jk) * vmask(ji  ,jj-1,jk) * e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk)
553               zet = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
554               zhdiv(ji,jj,jk) = ( zu - zu1 + zv - zv1 ) * zet 
555            END DO
556         END DO
557      END DO
558      CALL lbc_lnk( zhdiv, 'T', 1. )      ! Lateral boundary conditions on zhdiv
559      !
560      ! computation of vertical velocity from the bottom
561      pw(:,:,jpk) = 0._wp
562      DO jk = jpkm1, 1, -1
563         pw(:,:,jk) = pw(:,:,jk+1) - fse3t(:,:,jk) * zhdiv(:,:,jk)
564      END DO
565      !
566   END SUBROUTINE wzv
567
568
569   SUBROUTINE dta_eiv( kt )
570      !!----------------------------------------------------------------------
571      !!                  ***  ROUTINE dta_eiv  ***
572      !!
573      !! ** Purpose :   Compute the eddy induced velocity coefficient from the
574      !!      growth rate of baroclinic instability.
575      !!
576      !! ** Method : Specific to the offline model. Computes the horizontal
577      !!             values from the vertical value
578      !!----------------------------------------------------------------------
579      INTEGER, INTENT( in ) ::   kt     ! ocean time-step inedx
580      !!
581      INTEGER ::   ji, jj           ! dummy loop indices
582      !!----------------------------------------------------------------------
583      !
584      IF( kt == nit000 ) THEN
585         IF(lwp) WRITE(numout,*)
586         IF(lwp) WRITE(numout,*) 'dta_eiv : eddy induced velocity coefficients'
587         IF(lwp) WRITE(numout,*) '~~~~~~~'
588      ENDIF
589      !
590#if defined key_ldfeiv
591      ! Average the diffusive coefficient at u- v- points
592      DO jj = 2, jpjm1
593         DO ji = fs_2, fs_jpim1   ! vector opt.
594            aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) )
595            aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) )
596         END DO
597      END DO
598      CALL lbc_lnk( aeiu, 'U', 1. )   ;   CALL lbc_lnk( aeiv, 'V', 1. )    ! lateral boundary condition
599#endif
600      !
601   END SUBROUTINE dta_eiv
602
603
604   SUBROUTINE tau2wnd( ptaux, ptauy, pwspd )
605      !!---------------------------------------------------------------------
606      !!                    ***  ROUTINE sbc_tau2wnd  ***
607      !!
608      !! ** Purpose : Estimation of wind speed as a function of wind stress
609      !!
610      !! ** Method  : |tau|=rhoa*Cd*|U|^2
611      !!---------------------------------------------------------------------
612      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptaux, ptauy   ! wind stress in i-j direction resp.
613      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pwspd          ! wind speed
614      !!
615      REAL(wp) ::   zrhoa  = 1.22_wp       ! Air density kg/m3
616      REAL(wp) ::   zcdrag = 1.5e-3_wp     ! drag coefficient
617      REAL(wp) ::   ztx, zty, ztau, zcoef  ! temporary variables
618      INTEGER  ::   ji, jj                 ! dummy indices
619      !!---------------------------------------------------------------------
620      zcoef = 1. / ( zrhoa * zcdrag )
621!CDIR NOVERRCHK
622      DO jj = 2, jpjm1
623!CDIR NOVERRCHK
624         DO ji = fs_2, fs_jpim1   ! vector opt.
625            ztx = ptaux(ji,jj) * umask(ji,jj,1) + ptaux(ji-1,jj  ) * umask(ji-1,jj  ,1)
626            zty = ptauy(ji,jj) * vmask(ji,jj,1) + ptauy(ji  ,jj-1) * vmask(ji  ,jj-1,1)
627            ztau = 0.5 * SQRT( ztx * ztx + zty * zty )
628            pwspd(ji,jj) = SQRT ( ztau * zcoef ) * tmask(ji,jj,1)
629         END DO
630      END DO
631      CALL lbc_lnk( pwspd(:,:), 'T', 1. )
632      !
633   END SUBROUTINE tau2wnd
634
635
636   SUBROUTINE swap_dyn_data
637      !!----------------------------------------------------------------------
638      !!                    ***  ROUTINE swap_dyn_data  ***
639      !!
640      !! ** Purpose :   swap array data
641      !!----------------------------------------------------------------------
642      !
643      ! swap from record 2 to 1
644      tdta   (:,:,:,1) = tdta   (:,:,:,2)
645      sdta   (:,:,:,1) = sdta   (:,:,:,2)
646      avtdta (:,:,:,1) = avtdta (:,:,:,2)
647      udta   (:,:,:,1) = udta   (:,:,:,2)
648      vdta   (:,:,:,1) = vdta   (:,:,:,2)
649      wdta   (:,:,:,1) = wdta   (:,:,:,2)
650#if defined key_ldfslp && ! defined key_c1d
651      uslpdta (:,:,:,1) = uslpdta (:,:,:,2)
652      vslpdta (:,:,:,1) = vslpdta (:,:,:,2)
653      wslpidta(:,:,:,1) = wslpidta(:,:,:,2)
654      wslpjdta(:,:,:,1) = wslpjdta(:,:,:,2)
655#endif
656      hmlddta(:,:,1) = hmlddta(:,:,2) 
657      wspddta(:,:,1) = wspddta(:,:,2) 
658      frlddta(:,:,1) = frlddta(:,:,2) 
659      empdta (:,:,1) = empdta (:,:,2) 
660      qsrdta (:,:,1) = qsrdta (:,:,2) 
661      IF( l_offbbl ) THEN
662         bblxdta(:,:,1) = bblxdta(:,:,2)
663         bblydta(:,:,1) = bblydta(:,:,2) 
664      ENDIF
665
666#if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv
667      aeiwdta(:,:,1) = aeiwdta(:,:,2)
668#endif
669
670#if defined key_degrad
671      ahtudta(:,:,:,1) = ahtudta(:,:,:,2)
672      ahtvdta(:,:,:,1) = ahtvdta(:,:,:,2)
673      ahtwdta(:,:,:,1) = ahtwdta(:,:,:,2)
674#  if defined key_traldf_eiv
675      aeiudta(:,:,:,1) = aeiudta(:,:,:,2)
676      aeivdta(:,:,:,1) = aeivdta(:,:,:,2)
677      aeiwdta(:,:,:,1) = aeiwdta(:,:,:,2)
678#  endif
679#endif
680      !
681   END SUBROUTINE swap_dyn_data
682
683
684   SUBROUTINE assign_dyn_data
685      !!----------------------------------------------------------------------
686      !!                    ***  ROUTINE assign_dyn_data  ***
687      !!
688      !! ** Purpose :   Assign dynamical data to the data that have been read
689      !!                without time interpolation
690      !!
691      !!----------------------------------------------------------------------
692     
693      tsn(:,:,:,jp_tem) = tdta  (:,:,:,2)
694      tsn(:,:,:,jp_sal) = sdta  (:,:,:,2)
695      avt(:,:,:)        = avtdta(:,:,:,2)
696     
697      un (:,:,:) = udta  (:,:,:,2) 
698      vn (:,:,:) = vdta  (:,:,:,2)
699      wn (:,:,:) = wdta  (:,:,:,2)
700     
701#if defined key_ldfslp && ! defined key_c1d
702      uslp (:,:,:) = uslpdta (:,:,:,2) 
703      vslp (:,:,:) = vslpdta (:,:,:,2) 
704      wslpi(:,:,:) = wslpidta(:,:,:,2) 
705      wslpj(:,:,:) = wslpjdta(:,:,:,2) 
706#endif
707
708      hmld(:,:) = hmlddta(:,:,2) 
709      wndm(:,:) = wspddta(:,:,2) 
710      fr_i(:,:) = frlddta(:,:,2) 
711      emp (:,:) = empdta (:,:,2) 
712      emps(:,:) = emp(:,:) 
713      qsr (:,:) = qsrdta (:,:,2) 
714#if defined key_trabbl
715      IF( l_offbbl ) THEN
716         ahu_bbl(:,:) = bblxdta(:,:,2)
717         ahv_bbl(:,:) = bblydta(:,:,2) 
718      ENDIF
719#endif
720#if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv
721      aeiw(:,:) = aeiwdta(:,:,2)
722#endif
723     
724#if defined key_degrad
725      ahtu(:,:,:) = ahtudta(:,:,:,2)
726      ahtv(:,:,:) = ahtvdta(:,:,:,2)
727      ahtw(:,:,:) = ahtwdta(:,:,:,2)
728#  if defined key_traldf_eiv
729      aeiu(:,:,:) = aeiudta(:,:,:,2)
730      aeiv(:,:,:) = aeivdta(:,:,:,2)
731      aeiw(:,:,:) = aeiwdta(:,:,:,2)
732#  endif
733#endif
734      !
735   END SUBROUTINE assign_dyn_data
736
737
738   SUBROUTINE linear_interp_dyn_data( pweigh )
739      !!----------------------------------------------------------------------
740      !!               ***  ROUTINE linear_interp_dyn_data  ***
741      !!
742      !! ** Purpose :   linear interpolation of data
743      !!----------------------------------------------------------------------
744      REAL(wp), INTENT(in) ::   pweigh   ! weigh
745      !!
746      REAL(wp) :: zweighm1
747      !!----------------------------------------------------------------------
748
749      zweighm1 = 1. - pweigh
750     
751      tsn(:,:,:,jp_tem) = zweighm1 * tdta  (:,:,:,1) + pweigh * tdta  (:,:,:,2)
752      tsn(:,:,:,jp_sal) = zweighm1 * sdta  (:,:,:,1) + pweigh * sdta  (:,:,:,2)
753      avt(:,:,:)        = zweighm1 * avtdta(:,:,:,1) + pweigh * avtdta(:,:,:,2)
754     
755      un (:,:,:) = zweighm1 * udta  (:,:,:,1) + pweigh * udta  (:,:,:,2) 
756      vn (:,:,:) = zweighm1 * vdta  (:,:,:,1) + pweigh * vdta  (:,:,:,2)
757      wn (:,:,:) = zweighm1 * wdta  (:,:,:,1) + pweigh * wdta  (:,:,:,2)
758     
759#if defined key_ldfslp && ! defined key_c1d
760      uslp (:,:,:) = zweighm1 * uslpdta (:,:,:,1) + pweigh * uslpdta (:,:,:,2) 
761      vslp (:,:,:) = zweighm1 * vslpdta (:,:,:,1) + pweigh * vslpdta (:,:,:,2) 
762      wslpi(:,:,:) = zweighm1 * wslpidta(:,:,:,1) + pweigh * wslpidta(:,:,:,2) 
763      wslpj(:,:,:) = zweighm1 * wslpjdta(:,:,:,1) + pweigh * wslpjdta(:,:,:,2) 
764#endif
765
766      hmld(:,:) = zweighm1 * hmlddta(:,:,1) + pweigh  * hmlddta(:,:,2) 
767      wndm(:,:) = zweighm1 * wspddta(:,:,1) + pweigh  * wspddta(:,:,2) 
768      fr_i(:,:) = zweighm1 * frlddta(:,:,1) + pweigh  * frlddta(:,:,2) 
769      emp (:,:) = zweighm1 * empdta (:,:,1) + pweigh  * empdta (:,:,2) 
770      emps(:,:) = emp(:,:) 
771      qsr (:,:) = zweighm1 * qsrdta (:,:,1) + pweigh  * qsrdta (:,:,2) 
772#if defined key_trabbl
773      IF( l_offbbl ) THEN
774         ahu_bbl(:,:) = zweighm1 * bblxdta(:,:,1) +  pweigh  * bblxdta(:,:,2)
775         ahv_bbl(:,:) = zweighm1 * bblydta(:,:,1) +  pweigh  * bblydta(:,:,2)
776      ENDIF
777#endif
778
779#if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv
780      aeiw(:,:) = zweighm1 * aeiwdta(:,:,1) + pweigh * aeiwdta(:,:,2)
781#endif
782     
783#if defined key_degrad
784      ahtu(:,:,:) = zweighm1 * ahtudta(:,:,:,1) + pweigh * ahtudta(:,:,:,2)
785      ahtv(:,:,:) = zweighm1 * ahtvdta(:,:,:,1) + pweigh * ahtvdta(:,:,:,2)
786      ahtw(:,:,:) = zweighm1 * ahtwdta(:,:,:,1) + pweigh * ahtwdta(:,:,:,2)
787#  if defined key_traldf_eiv
788      aeiu(:,:,:) = zweighm1 * aeiudta(:,:,:,1) + pweigh * aeiudta(:,:,:,2)
789      aeiv(:,:,:) = zweighm1 * aeivdta(:,:,:,1) + pweigh * aeivdta(:,:,:,2)
790      aeiw(:,:,:) = zweighm1 * aeiwdta(:,:,:,1) + pweigh * aeiwdta(:,:,:,2)
791#  endif
792#endif
793      !     
794   END SUBROUTINE linear_interp_dyn_data
795
796   !!======================================================================
797END MODULE dtadyn
Note: See TracBrowser for help on using the repository browser.