source: trunk/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 10 years ago

Update NEMOGCM from branch nemo_v3_3_beta

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