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/nemo_v3_3_beta/NEMOGCM/NEMO/OFF_SRC – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90 @ 2445

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

suppression of line commenst

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