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 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 14 years ago

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