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.
hltinc.F90 in branches/TAM_V3_0/NEMOTAM/OPATAM_SRC – NEMO

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/hltinc.F90 @ 2587

Last change on this file since 2587 was 2587, checked in by vidard, 13 years ago

refer to ticket #798

File size: 26.4 KB
Line 
1MODULE hltinc
2#if defined key_tam
3   !!======================================================================
4   !!                       ***  MODULE hltinc  ***
5   !! Linear-Tangent test increment : Apply an increment to the control vect
6   !!                                 Increment is either generated by 
7   !!                                 difference between two states or either
8   !!                                 directly read from an input file
9   !!======================================================================
10
11   !!----------------------------------------------------------------------
12   !!   'key_asminc' : Switch on the assimilation increment interface
13   !!----------------------------------------------------------------------
14   !!   asm_inc_init : Initialize the increment arrays and IAU weights
15   !!   calc_date    : Compute the calendar date YYYYMMDD on a given step
16   !!   tra_asm_inc  : Apply the tracer (T and S) increments
17   !!   dyn_asm_inc  : Apply the dynamic (u and v) increments
18   !!   ssh_asm_inc  : Apply the SSH increment
19   !!----------------------------------------------------------------------
20   !! * Modules used   
21   USE par_kind, ONLY : &       ! Precision variables
22      & wp
23   USE in_out_manager, ONLY : & ! I/O manager
24      & lwp,      &
25      & numnam,   &
26      & numout,   &
27      & ctl_warn, &
28      & ctl_stop, &
29      & nit000,   &
30      & nstop,    &
31      & ln_rstart
32   USE par_oce, ONLY : & ! Ocean space and time domain variables
33      & jpi,  & 
34      & jpj,  &
35      & jpk,  &
36      & jpkm1
37   USE dom_oce, ONLY : & ! Ocean space and time domain
38      & rdt,    &
39      & n_cla,  &
40      & neuler, &
41      & ln_zps, &
42      & tmask,  &
43      & umask,  &
44      & vmask,  &
45      & nldi,   &
46      & nldj,   &
47      & nlei,   &
48      & nlej
49   USE c1d, ONLY : &    ! 1D initialization
50      & lk_c1d
51   USE oce, ONLY : &      ! Dynamics and active tracers defined in memory
52      & ub, un, ua,    &
53      & vb, vn, va,    &
54      & tb, tn, ta,    &
55      & sb, sn, sa,    &
56      & sshb, sshn,    &
57      & rhd, rhop,     &
58      & rotb, rotn,    &
59      & hdivb, hdivn,  &
60      & gtu, gsu, gru, &
61      & gtv, gsv, grv 
62   USE oce_tam, ONLY : &      ! Dynamics and active tracers defined in memory
63      & ub_tl, un_tl, ua_tl,    &
64      & vb_tl, vn_tl, va_tl,    &
65      & tb_tl, tn_tl, ta_tl,    &
66      & sb_tl, sn_tl, sa_tl,    &
67      & sshb_tl, sshn_tl,    &
68      & rhd_tl, rhop_tl,     &
69      & rotb_tl, rotn_tl,    &
70      & hdivb_tl, hdivn_tl,  &
71      & gtu_tl, gsu_tl, gru_tl, &
72      & gtv_tl, gsv_tl, grv_tl
73   USE divcur, ONLY : &   ! Horizontal divergence and relative vorticity
74      & div_cur
75   USE cla_div, ONLY : &  ! Specific update of the horizontal divergence
76      & div_cla           ! (specific to ORCA_R2)
77   USE wzvmod, ONLY : &   ! Vertical velocity
78      & wzv
79   USE eosbn2, ONLY : &   ! Equation of state - in situ and potential density
80      & eos
81   USE zpshde, ONLY : &   ! Partial step : Horizontal Derivative
82      & zps_hde 
83   USE divcur_tam, ONLY : &   ! Horizontal divergence and relative vorticity
84      & div_cur_tan
85   USE cla_div_tam, ONLY : &  ! Specific update of the horizontal divergence
86      & div_cla_tan           ! (specific to ORCA_R2)
87   USE wzvmod_tam, ONLY : &   ! Vertical velocity
88      & wzv_tan
89   USE eosbn2_tam, ONLY : &   ! Equation of state - in situ and potential density
90      & eos_tan
91   USE zpshde_tam, ONLY : &   ! Partial step : Horizontal Derivative
92      & zps_hde_tan
93   USE iom                ! Library to read input files
94   USE par_tlm
95
96   IMPLICIT NONE
97
98   !! * Routine accessibility
99   PRIVATE
100   PUBLIC hlt_inc_bld      !: Initialize the increment arrays
101
102   !! * Private Module variables
103   REAL(wp), PRIVATE, DIMENSION(:,:,:), ALLOCATABLE :: &
104      & t_hltinc1, &       !: Increment to the background temperature
105      & s_hltinc1, &       !: Increment to the background salinity
106      & u_hltinc1, &       !: Increment to the u-component velocity
107      & v_hltinc1, &       !: Increment to the v-component velocity
108      & t_hltinc2, &       !: Increment to the background temperature
109      & s_hltinc2, &       !: Increment to the background salinity
110      & u_hltinc2, &       !: Increment to the u-component velocity
111      & v_hltinc2          !: Increment to the v-component velocity
112
113#if defined key_dynspg_flt         
114   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: &
115      & ssh_hltinc1, &     !: Increment to the background sea surface height
116      & ssh_hltinc2        !: Increment to the background sea surface height
117#endif
118   CHARACTER (LEN=40), PRIVATE, PARAMETER :: &
119      & c_hltincwri = 'increments_01'                 !: Filename for storing the
120                                                   !: hlt increment dX
121
122CONTAINS
123
124   SUBROUTINE hlt_inc_bld ( pstg )
125      !!----------------------------------------------------------------------
126      !!                    ***  ROUTINE hlt_inc_bld  ***
127      !!         
128      !! ** Purpose : Initialize the assimilation increment and IAU weights.
129      !!
130      !! ** Method  : Initialize the assimilation increment and IAU weights.
131      !!
132      !! ** Action  :
133      !!
134      !! History :
135      !!        !  10-07  (F. Vigilant) Original code
136      !!----------------------------------------------------------------------
137
138      IMPLICIT NONE
139      INTEGER, INTENT(IN) :: &
140         & pstg               ! Current stage
141      !! * Modules used
142      !! * Local declarations
143      INTEGER :: &
144         & jk, jj, ji
145      INTEGER :: &
146         & inum
147      REAL(wp), DIMENSION(5) :: &
148         & zvalmax
149      REAL(wp) :: &
150         & zt, zs, zu, zv, zssh, &
151         & znormt, znorms, znormu, znormv, znormssh
152      REAL(wp) :: &
153         & znorm
154      LOGICAL :: &
155         & zlnrm = .FALSE.
156
157      !--------------------------------------------------------------------
158      ! Initialize the Incremental Analysis Updating weighting function
159      !--------------------------------------------------------------------
160
161      IF ( .NOT. ln_incdx ) THEN      ! No increment file defined
162
163         c_hltrst1 = 'restart.nc'
164         c_hltrst2 = 'restart2.nc'
165
166          IF(lwp) THEN
167             WRITE(numout,*)
168             WRITE(numout,*) 'hlt_inc_init : Open restart file 1 ', c_hltrst1
169             WRITE(numout,*) '~~~~~~~~~~~~'
170          ENDIF
171
172          !--------------------------------------------------------------------
173          ! Allocate and initialize the increment arrays 1
174          !--------------------------------------------------------------------
175          ALLOCATE( t_hltinc1(jpi,jpj,jpk) )
176          ALLOCATE( s_hltinc1(jpi,jpj,jpk) )
177          ALLOCATE( u_hltinc1(jpi,jpj,jpk) )
178          ALLOCATE( v_hltinc1(jpi,jpj,jpk) )
179#if defined key_dynspg_flt
180          ALLOCATE( ssh_hltinc1(jpi,jpj  ) )
181#endif
182          t_hltinc1(:,:,:) = 0.0_wp
183          s_hltinc1(:,:,:) = 0.0_wp
184          u_hltinc1(:,:,:) = 0.0_wp
185          v_hltinc1(:,:,:) = 0.0_wp
186#if defined key_dynspg_flt
187          ssh_hltinc1(:,:) = 0.0_wp
188#endif
189         CALL iom_open( c_hltrst1, inum )
190
191         IF ( ln_hltt  )   CALL iom_get( inum, jpdom_autoglo, 'tn',   t_hltinc1 ) 
192         IF ( ln_hlts  )   CALL iom_get( inum, jpdom_autoglo, 'sn',   s_hltinc1 )
193         IF ( ln_hltuv )   CALL iom_get( inum, jpdom_autoglo, 'un',   u_hltinc1 )
194         IF ( ln_hltuv )   CALL iom_get( inum, jpdom_autoglo, 'vn',   v_hltinc1 )
195#if defined key_dynspg_flt
196         IF ( ln_hltuv )   CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_hltinc1 )
197#endif
198         CALL iom_close( inum )
199
200          IF(lwp) THEN
201             WRITE(numout,*)
202             WRITE(numout,*) 'hlt_inc_init : Open restart file 2 ', c_hltrst2
203             WRITE(numout,*) '~~~~~~~~~~~~'
204          ENDIF
205
206          !--------------------------------------------------------------------
207          ! Allocate and initialize the increment arrays 2
208          !--------------------------------------------------------------------
209          ALLOCATE( t_hltinc2(jpi,jpj,jpk) )
210          ALLOCATE( s_hltinc2(jpi,jpj,jpk) )
211          ALLOCATE( u_hltinc2(jpi,jpj,jpk) )
212          ALLOCATE( v_hltinc2(jpi,jpj,jpk) )
213#if defined key_dynspg_flt
214          ALLOCATE( ssh_hltinc2(jpi,jpj  ) )
215#endif
216          t_hltinc2(:,:,:) = 0.0_wp
217          s_hltinc2(:,:,:) = 0.0_wp
218          u_hltinc2(:,:,:) = 0.0_wp
219          v_hltinc2(:,:,:) = 0.0_wp
220#if defined key_dynspg_flt
221          ssh_hltinc2(:,:) = 0.0_wp
222#endif
223         CALL iom_open( c_hltrst2, inum )
224
225         IF ( ln_hltt  )   CALL iom_get( inum, jpdom_autoglo, 'tn',   t_hltinc2 ) 
226         IF ( ln_hlts  )   CALL iom_get( inum, jpdom_autoglo, 'sn',   s_hltinc2 )
227         IF ( ln_hltuv )   CALL iom_get( inum, jpdom_autoglo, 'un',   u_hltinc2 )
228         IF ( ln_hltuv )   CALL iom_get( inum, jpdom_autoglo, 'vn',   v_hltinc2 )
229#if defined key_dynspg_flt
230         IF ( ln_hltuv )   CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_hltinc2 )
231#endif
232         CALL iom_close( inum )
233
234          !--------------------------------------------------------------------
235          ! Apply difference restart1 - restart2
236          !--------------------------------------------------------------------
237          t_hltinc1(:,:,:) = t_hltinc1(:,:,:) - t_hltinc2(:,:,:)
238          s_hltinc1(:,:,:) = s_hltinc1(:,:,:) - s_hltinc2(:,:,:)
239          u_hltinc1(:,:,:) = u_hltinc1(:,:,:) - u_hltinc2(:,:,:)
240          v_hltinc1(:,:,:) = v_hltinc1(:,:,:) - v_hltinc2(:,:,:)
241#if defined key_dynspg_flt
242          ssh_hltinc1(:,:) = ssh_hltinc1(:,:) - ssh_hltinc2(:,:)
243#endif
244          !--------------------------------------------------------------------
245          ! Deallocate and initialize the increment arrays 2
246          !--------------------------------------------------------------------
247          DEALLOCATE( t_hltinc2 )
248          DEALLOCATE( s_hltinc2 )
249          DEALLOCATE( u_hltinc2 )
250          DEALLOCATE( v_hltinc2 )
251#if defined key_dynspg_flt
252          DEALLOCATE( ssh_hltinc2 )
253#endif
254      ELSE    ! increment file defined
255
256!          c_hltinc = c_hltincwri
257          WRITE(c_hltinc, FMT='(A,A)' ) TRIM( c_hltincwri ), '.nc'
258          IF(lwp) THEN
259             WRITE(numout,*)
260             WRITE(numout,*) 'hlt_inc_init : Open increment file ', c_hltinc
261             WRITE(numout,*) '~~~~~~~~~~~~'
262          ENDIF
263
264          !--------------------------------------------------------------------
265          ! Allocate and initialize the increment arrays 1
266          !--------------------------------------------------------------------
267          ALLOCATE( t_hltinc1(jpi,jpj,jpk) )
268          ALLOCATE( s_hltinc1(jpi,jpj,jpk) )
269          ALLOCATE( u_hltinc1(jpi,jpj,jpk) )
270          ALLOCATE( v_hltinc1(jpi,jpj,jpk) )
271#if defined key_dynspg_flt
272          ALLOCATE( ssh_hltinc1(jpi,jpj)   )
273#endif
274          t_hltinc1(:,:,:) = 0.0_wp
275          s_hltinc1(:,:,:) = 0.0_wp
276          u_hltinc1(:,:,:) = 0.0_wp
277          v_hltinc1(:,:,:) = 0.0_wp
278          ssh_hltinc1(:,:) = 0.0_wp
279
280         CALL iom_open( c_hltinc, inum )
281
282         IF ( ln_hltt   ) CALL iom_get( inum, jpdom_autoglo, 'bckint',   t_hltinc1 ) 
283         IF ( ln_hlts   ) CALL iom_get( inum, jpdom_autoglo, 'bckins',   s_hltinc1 )
284         IF ( ln_hltuv  ) CALL iom_get( inum, jpdom_autoglo, 'bckinu',   u_hltinc1 )
285         IF ( ln_hltuv  ) CALL iom_get( inum, jpdom_autoglo, 'bckinv',   v_hltinc1 )
286#if defined key_dynspg_flt
287         IF ( ln_hltssh ) CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_hltinc1 )
288#endif
289         CALL iom_close( inum )
290         
291      ENDIF
292
293      !--------------------------------------------------------------------
294      ! Check X max value and normalize if needed
295      !--------------------------------------------------------------------
296
297      IF ( ln_hnorm ) THEN
298
299         zvalmax(1) = maxval(ABS(t_hltinc1(:,:,:)))
300         zvalmax(2) = maxval(ABS(s_hltinc1(:,:,:)))
301         zvalmax(3) = maxval(ABS(u_hltinc1(:,:,:)))
302         zvalmax(4) = maxval(ABS(v_hltinc1(:,:,:)))
303#if defined key_dynspg_flt
304         IF ( ln_incdx ) THEN 
305            zvalmax(5) = maxval(ABS(ssh_hltinc1(:,:)))
306         ENDIF
307#endif
308         IF(lwp) THEN
309            WRITE(numout,*)
310            WRITE(numout,*) '        Increments max value before normalization'
311            WRITE(numout,*) '        Temperature increment max value         : ', zvalmax(1)
312            WRITE(numout,*) '        Salinity increment max value            : ', zvalmax(2)
313            WRITE(numout,*) '        Horizontal velocity increment max value : ', zvalmax(3)
314            WRITE(numout,*) '        Zonal velocity increment max value      : ', zvalmax(4)
315            IF ( ln_incdx ) THEN
316               WRITE(numout,*) '        Sea Surface increment max value         : ', zvalmax(5)
317            ENDIF
318            WRITE(numout,*) 
319            WRITE(numout,*) '        Normalization bounds:'
320            WRITE(numout,*) '        Temperature rhstdt                      : ',rhstdt
321            WRITE(numout,*) '        Salinity    rhstds                      : ',rhstds
322            WRITE(numout,*) '        Velocity    rhstduv                     : ',rhstduv
323#if defined key_dynspg_flt
324            IF ( ln_incdx ) THEN
325               WRITE(numout,*) '        Sea Surface Height    rhstdssh          : ',rhstdssh
326            ENDIF
327#endif
328         ENDIF
329         
330         IF ( zvalmax(1) > rhstdt )   THEN
331            zlnrm = .TRUE.
332         ELSE
333            zvalmax(1) = rhstdt
334         ENDIF
335         IF ( zvalmax(2) > rhstds )   THEN
336            zlnrm = .TRUE.
337         ELSE
338            zvalmax(2) = rhstds
339         ENDIF
340         IF ( zvalmax(3) > rhstduv )   THEN
341            zlnrm = .TRUE.
342         ELSE
343            zvalmax(3) = rhstduv
344         ENDIF
345         IF ( zvalmax(4) > rhstduv )   THEN
346            zlnrm = .TRUE.
347         ELSE
348            zvalmax(4) = rhstduv
349         ENDIF
350#if defined key_dynspg_flt
351         IF ( ln_incdx ) THEN
352            IF ( zvalmax(5) > rhstdssh )   THEN
353               zlnrm = .TRUE.
354            ELSE
355               zvalmax(5) = rhstdssh
356            ENDIF
357         ENDIF
358#endif
359         
360         IF ( zlnrm ) THEN
361!!$         IF ( .NOT. ln_incdx ) THEN
362!!$            znorm = MIN( stdt/zvalmax(1), stds/zvalmax(2), stdu/zvalmax(3), &
363!!$                       & stdv/zvalmax(4), stdssh/zvalmax(5))
364!!$         ELSE
365!!$            znorm = MIN( stdt/zvalmax(1), stds/zvalmax(2), stdu/zvalmax(3), &
366!!$                       & stdv/zvalmax(4))
367!!$
368!!$         ENDIF
369            znormt = rhstdt/zvalmax(1)
370            znorms = rhstds/zvalmax(2)
371            znormu = rhstduv/zvalmax(3)
372            znormv = rhstduv/zvalmax(4)
373#if defined key_dynspg_flt
374            IF ( ln_incdx ) znormssh = rhstdssh/zvalmax(5)
375#endif           
376            DO jk = 1, jpk
377               DO jj = nldj, nlej
378                  DO ji = nldi, nlei
379                     t_hltinc1(ji,jj,jk) = t_hltinc1(ji,jj,jk) * znormt
380                  END DO
381               END DO
382            END DO
383            DO jk = 1, jpk
384               DO jj = nldj, nlej
385                  DO ji = nldi, nlei
386                     s_hltinc1(ji,jj,jk) = s_hltinc1(ji,jj,jk) * znorms
387                  END DO
388               END DO
389            END DO
390            DO jk = 1, jpk
391               DO jj = nldj, nlej
392                  DO ji = nldi, nlei
393                     u_hltinc1(ji,jj,jk) = u_hltinc1(ji,jj,jk) * znormu
394                  END DO
395               END DO
396            END DO
397            DO jk = 1, jpk
398               DO jj = nldj, nlej
399                  DO ji = nldi, nlei
400                     v_hltinc1(ji,jj,jk) = v_hltinc1(ji,jj,jk) * znormv
401                  END DO
402               END DO
403            END DO
404#if defined key_dynspg_flt
405            IF ( ln_incdx ) THEN
406               DO jj = nldj, nlej
407                  DO ji = nldi, nlei
408                     ssh_hltinc1(ji,jj) = ssh_hltinc1(ji,jj) * znormssh
409                  END DO
410               END DO
411            ENDIF
412#endif           
413         ENDIF
414
415      ENDIF
416
417      zt = maxval(ABS(t_hltinc1(:,:,:)))
418      zs = maxval(ABS(s_hltinc1(:,:,:)))
419      zu = maxval(ABS(u_hltinc1(:,:,:)))
420      zv = maxval(ABS(v_hltinc1(:,:,:)))
421#if defined key_dynspg_flt
422        zssh =  maxval(ABS(ssh_hltinc1(:,:)))
423#endif
424      IF(lwp) THEN
425         WRITE(numout,*)
426         WRITE(numout,*) '        Increments max value after normalization'
427         WRITE(numout,*) '        Temperature increment max value         : ', zt
428         WRITE(numout,*) '        Salinity increment max value            : ', zs
429         WRITE(numout,*) '        Horizontal velocity increment max value : ', zu
430         WRITE(numout,*) '        Zonal velocity increment max value      : ', zv
431#if defined key_dynspg_flt
432            WRITE(numout,*) '        Sea Surface increment max value         : ', zssh
433#endif
434         WRITE(numout,*) 
435      ENDIF
436
437      SELECT CASE( nstg )
438      CASE ( 1 )                                 ! X = X + dX
439         IF ( .NOT. ln_incdx ) CALL hlt_inc_wri( ln_incdx, nstg)
440                               CALL tra_hlt_inc            ! Tracers
441                               CALL dyn_hlt_inc            ! Dynamics
442#if defined key_dynspg_flt
443         IF ( ln_incdx )       CALL ssh_hlt_inc            ! SSH
444#endif
445      CASE ( 2 )                                 ! X_tl = dX
446         IF ( .NOT. ln_incdx ) CALL hlt_inc_wri( ln_incdx, nstg)
447                               CALL tra_hlt_inc_tan        ! Tracers
448                               CALL dyn_hlt_inc_tan        ! Dynamics
449#if defined key_dynspg_flt
450                               CALL ssh_hlt_inc_tan        ! SSH
451#endif
452      CASE ( 3 )            ! save dX to file if compute from restart
453         IF ( .NOT. ln_incdx ) CALL hlt_inc_wri( ln_incdx, nstg)
454      END SELECT
455
456   END SUBROUTINE hlt_inc_bld
457
458   SUBROUTINE tra_hlt_inc
459      !!----------------------------------------------------------------------
460      !!                    ***  ROUTINE tra_hlt_inc  ***
461      !!         
462      !! ** Purpose : Apply the tracer (T and S)  increments
463      !!
464      !! ** Method  :
465      !!
466      !! ** Action  :
467      !!
468      !! History :
469      !!        !  10-07  (F. Vigilant) Original code
470      !!----------------------------------------------------------------------
471
472      IMPLICIT NONE
473
474      !! * Arguments
475      !! * Local declarations
476           
477      neuler = 0  ! Force Euler forward step
478
479      ! Add increment
480      tn(:,:,:) = tn(:,:,:) + t_hltinc1(:,:,:)   
481      sn(:,:,:) = sn(:,:,:) + s_hltinc1(:,:,:)   
482
483      tb(:,:,:) = tn(:,:,:)                        ! Update before fields
484      sb(:,:,:) = sn(:,:,:)
485
486      CALL eos( tb, sb, rhd, rhop )                ! Before potential and in situ densities
487
488      IF( ln_zps .AND. .NOT. lk_c1d ) &
489          &  CALL zps_hde( nit000, tb, sb, rhd,  &  ! Partial steps: before horizontal derivative
490          &                gtu, gsu, gru,        &  ! of T, S, rd at the bottom ocean level
491          &                gtv, gsv, grv )
492
493      DEALLOCATE( t_hltinc1 )
494      DEALLOCATE( s_hltinc1 )
495
496   END SUBROUTINE tra_hlt_inc
497
498   SUBROUTINE dyn_hlt_inc
499      !!----------------------------------------------------------------------
500      !!                    ***  ROUTINE dyn_hlt_inc  ***
501      !!         
502      !! ** Purpose : Apply the dynamics (u and v) increments.
503      !!
504      !! ** Method  :
505      !!
506      !! ** Action  :
507      !!
508      !! History :
509      !!        !  10-07  (F. Vigilant) Original code
510      !!----------------------------------------------------------------------
511
512      IMPLICIT NONE
513
514      !! * Arguments
515      !! * Local declarations
516      INTEGER :: &
517         & kt
518
519      kt = nit000
520      neuler = 0                    ! Force Euler forward step
521
522      ! Add increment
523      un(:,:,:) = un(:,:,:) + u_hltinc1(:,:,:)
524      vn(:,:,:) = vn(:,:,:) + v_hltinc1(:,:,:) 
525
526      ub(:,:,:) = un(:,:,:)         ! Update before fields
527      vb(:,:,:) = vn(:,:,:)
528 
529      CALL div_cur( kt )                  ! Compute divergence and curl for now fields
530      IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence)
531
532      rotb (:,:,:) = rotn (:,:,:)   ! Update before fields
533      hdivb(:,:,:) = hdivn(:,:,:)
534
535      CALL wzv( kt )                ! Vertical velocity
536
537      DEALLOCATE( u_hltinc1 )
538      DEALLOCATE( v_hltinc1 )
539
540   END SUBROUTINE dyn_hlt_inc
541
542   SUBROUTINE ssh_hlt_inc
543      !!----------------------------------------------------------------------
544      !!                    ***  ROUTINE ssh_hlt_inc  ***
545      !!         
546      !! ** Purpose : Apply the sea surface height increment.
547      !!
548      !! ** Method  :
549      !!
550      !! ** Action  :
551      !!
552      !! History :
553      !!        !  10-07  (F. Vigilant) Original code
554      !!----------------------------------------------------------------------
555
556      IMPLICIT NONE
557
558      !! * Arguments
559
560      !! * Local declarations
561      INTEGER :: &
562         & it
563
564      neuler = 0                    ! Force Euler forward step
565
566      ! Add increment
567      sshn(:,:) = sshn(:,:) + ssh_hltinc1(:,:) 
568
569      sshb(:,:) = sshn(:,:)         ! Update before fields
570
571      DEALLOCATE( ssh_hltinc1 )
572
573   END SUBROUTINE ssh_hlt_inc
574
575   SUBROUTINE tra_hlt_inc_tan
576      !!----------------------------------------------------------------------
577      !!                    ***  ROUTINE tra_hlt_inc_tan  ***
578      !!         
579      !! ** Purpose : Apply the tracer (T and S)  increments
580      !!
581      !! ** Method  :
582      !!
583      !! ** Action  :
584      !!
585      !! History :
586      !!        !  10-07  (F. Vigilant) Original code
587      !!----------------------------------------------------------------------
588
589      IMPLICIT NONE
590
591      !! * Arguments
592      !! * Local declarations
593           
594      neuler = 0  ! Force Euler forward step
595
596      ! Add increment
597      tn_tl(:,:,:) = t_hltinc1(:,:,:)   
598      sn_tl(:,:,:) = s_hltinc1(:,:,:)   
599
600      CALL lbc_lnk( tn_tl, 'T',  1.0 )
601      CALL lbc_lnk( sn_tl, 'T',  1.0 )
602
603      tb_tl(:,:,:) = tn_tl(:,:,:)                        ! Update before fields
604      sb_tl(:,:,:) = sn_tl(:,:,:)
605
606      CALL eos_tan( tb, sb, tb_tl, sb_tl, rhd_tl, rhop_tl )                ! Before potential and in situ densities
607         
608      IF( ln_zps .AND. .NOT. lk_c1d ) &
609          &  CALL zps_hde_tan( nit000, tb, sb, tb_tl, sb_tl, rhd_tl,  &  ! Partial steps: before horizontal derivative
610          &                gtu_tl, gsu_tl, gru_tl,        &  ! of T, S, rd at the bottom ocean level
611          &                gtv_tl, gsv_tl, grv_tl )
612
613      DEALLOCATE( t_hltinc1 )
614      DEALLOCATE( s_hltinc1 )
615
616   END SUBROUTINE tra_hlt_inc_tan
617
618   SUBROUTINE dyn_hlt_inc_tan
619      !!----------------------------------------------------------------------
620      !!                    ***  ROUTINE dyn_hlt_inc_tan  ***
621      !!         
622      !! ** Purpose : Apply the dynamics (u and v) increments.
623      !!
624      !! ** Method  :
625      !!
626      !! ** Action  :
627      !!
628      !! History :
629      !!        !  10-07  (F. Vigilant) Original code
630      !!----------------------------------------------------------------------
631
632      IMPLICIT NONE
633
634      !! * Arguments
635      !! * Local declarations
636      INTEGER :: &
637         & kt
638
639      kt = nit000
640      neuler = 0                    ! Force Euler forward step
641
642      ! Add increment
643      un_tl(:,:,:) = u_hltinc1(:,:,:)
644      vn_tl(:,:,:) = v_hltinc1(:,:,:) 
645
646      CALL lbc_lnk( un_tl, 'U', -1.0 )
647      CALL lbc_lnk( vn_tl, 'V', -1.0 )
648 
649      ub_tl(:,:,:) = un_tl(:,:,:)         ! Update before fields
650      vb_tl(:,:,:) = vn_tl(:,:,:)
651 
652      CALL div_cur_tan( kt )                  ! Compute divergence and curl for now fields
653      IF( n_cla == 1 ) CALL div_cla_tan( kt ) ! Cross Land Advection (Update Hor. divergence)
654
655      rotb_tl (:,:,:) = rotn_tl (:,:,:)   ! Update before fields
656      hdivb_tl(:,:,:) = hdivn_tl(:,:,:)
657
658      CALL wzv_tan( kt )                ! Vertical velocity
659
660      DEALLOCATE( u_hltinc1 )
661      DEALLOCATE( v_hltinc1 )
662
663   END SUBROUTINE dyn_hlt_inc_tan
664
665   SUBROUTINE ssh_hlt_inc_tan
666      !!----------------------------------------------------------------------
667      !!                    ***  ROUTINE ssh_hlt_inc_tan  ***
668      !!         
669      !! ** Purpose : Apply the sea surface height increment.
670      !!
671      !! ** Method  :
672      !!
673      !! ** Action  :
674      !!
675      !! History :
676      !!        !  10-07  (F. Vigilant) Original code
677      !!----------------------------------------------------------------------
678
679      IMPLICIT NONE
680
681      !! * Arguments
682
683      !! * Local declarations
684      INTEGER :: &
685         & it
686
687      neuler = 0                    ! Force Euler forward step
688
689      ! Add increment
690      sshn_tl(:,:) = ssh_hltinc1(:,:) 
691
692      sshb_tl(:,:) = sshn_tl(:,:)         ! Update before fields
693
694      DEALLOCATE( ssh_hltinc1 )
695
696   END SUBROUTINE ssh_hlt_inc_tan
697
698   SUBROUTINE hlt_inc_wri ( plrst, pstg ) 
699      !!----------------------------------------------------------------------
700      !!                    ***  ROUTINE hlt_inc_wri  ***
701      !!         
702      !! ** Purpose : Write the  increment.
703      !!
704      !! ** Method  :
705      !!
706      !! ** Action  :
707      !!
708      !! History :
709      !!        !  10-07  (F. Vigilant) Original code
710      !!----------------------------------------------------------------------
711
712      !! *Module udes
713      USE iom 
714      !! * Arguments
715      LOGICAL, INTENT(in) :: &
716         & plrst            ! logical to write or not ssh
717      !! * Local declarations
718      INTEGER :: &
719         & inum, &                  ! File unit number
720         & fd                       ! field number
721      INTEGER :: &
722         & it, pstg
723      CHARACTER (LEN=100) :: &
724         & cl_hltinc         
725
726      cl_hltinc = c_hltincwri
727      fd = 1
728
729      WRITE(cl_hltinc, FMT='(A,A)' ) TRIM( cl_hltinc ), '_output.nc'
730!      WRITE(cl_hltinc, FMT='(A)' ) TRIM( cl_hltinc )
731      cl_hltinc = TRIM( cl_hltinc )
732      CALL iom_open( cl_hltinc, inum, ldwrt = .TRUE., kiolib = jprstlib)
733
734      IF(lwp) THEN
735         WRITE(numout,*)
736         WRITE(numout,*) 'Writing increment in file : ', c_hltinc
737         WRITE(numout,*)
738      ENDIF
739
740      ! Output increment fields
741      CALL iom_rstput( fd, fd, inum, 'bckinu'   , u_hltinc1   )
742      CALL iom_rstput( fd, fd, inum, 'bkcinv'   , v_hltinc1   )
743      CALL iom_rstput( fd, fd, inum, 'bckint'   , t_hltinc1   )
744      CALL iom_rstput( fd, fd, inum, 'bckins'   , s_hltinc1   )
745#if defined key_dynspg_flt
746      IF ( plrst ) THEN
747         CALL iom_rstput( fd, fd, inum, 'bckineta' , ssh_hltinc1 )
748      ENDIF
749#endif
750      CALL iom_close( inum )
751
752      IF ( pstg == 3 ) THEN
753         DEALLOCATE( t_hltinc1 )
754         DEALLOCATE( s_hltinc1 )
755         DEALLOCATE( u_hltinc1 )
756         DEALLOCATE( v_hltinc1 )
757#if defined key_dynspg_flt
758         IF (  plrst ) THEN
759            DEALLOCATE( ssh_hltinc1 )
760         ENDIF
761#endif
762      ENDIF
763
764   END SUBROUTINE hlt_inc_wri
765#endif
766END MODULE hltinc
Note: See TracBrowser for help on using the repository browser.