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.
sbcfwb_tam.F90 in branches/2010_and_older/TAM_V3_2_2/NEMOTAM/OPATAM_SRC/SBC – NEMO

source: branches/2010_and_older/TAM_V3_2_2/NEMOTAM/OPATAM_SRC/SBC/sbcfwb_tam.F90 @ 7797

Last change on this file since 7797 was 2578, checked in by rblod, 13 years ago

first import of NEMOTAM 3.2.2

File size: 18.9 KB
Line 
1MODULE sbcfwb_tam
2#if defined key_tam
3   !!======================================================================
4   !!                       ***  MODULE  sbcfwb  ***
5   !! Ocean fluxes   : domain averaged freshwater budget
6   !!======================================================================
7   !! History of the direct module:
8   !!            8.2  !  01-02  (E. Durand)  Original code
9   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
10   !!            9.0  !  06-08  (G. Madec)  Surface module
11   !!            9.2  !  09-07  (C. Talandier) emp mean s spread over erp area
12   !! History of the T&A module:
13   !!            9.0  !  06-08  (A. Vidard)  Surface module
14   !!            9.2  !  09-07  (A. Vidard)  NEMO3.2 update
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   sbc_fwb      : freshwater budget for global ocean configurations
19   !!                  in free surface and forced mode
20   !!----------------------------------------------------------------------
21   USE par_kind      , ONLY: & !
22      & wp
23   USE par_oce       , ONLY: & ! Ocean space and time domain variables
24      & jpi,                 &
25      & jpj,                 &
26      & jpiglo
27   USE oce_tam       , ONLY: & ! ocean dynamics and tracers
28      & sshn_tl,             &
29      & sshn_ad
30   USE sbc_oce_tam   , ONLY: & ! surface variables
31      & emp_tl,              &
32      & emps_tl,             &
33      & emp_ad,              &
34      & emps_ad
35   USE dom_oce       , ONLY: & ! ocean space and time domain
36      & rdttra,              &
37      & e1t,                 &
38      & e2t,                 &
39#if defined key_zco
40      & e3t_0,               &
41#else
42      & e3t,                 &
43#endif
44      & tmask,               &
45      & tmask_i,             &
46      & mig,                 &
47      & mjg,                 &
48      & nldi,                &
49      & nldj,                &
50      & nlei,                &
51      & nlej
52
53   USE tstool_tam    , ONLY: &
54      & prntst_adj,          &
55      & stdemp,              & !   evaporation minus precip
56      & stdssh                 !   sea surface height
57   USE in_out_manager, ONLY: & ! I/O manager
58      & lwp,                 &
59      & ctl_stop,            &
60      & numout,              & 
61      & nit000,              & 
62      & nitend,              & 
63      & ctmp1 
64   USE lib_mpp                 ! distribued memory computing library
65   USE gridrandom    , ONLY: & ! Random Gaussian noise on grids
66      & grid_random
67   USE dotprodfld    , ONLY: & ! Computes dot product for 3D and 2D fields
68      & dot_product
69
70   IMPLICIT NONE
71   PRIVATE
72
73   PUBLIC   sbc_fwb_adj      ! routine called by sbcmod_tam
74   PUBLIC   sbc_fwb_tan      ! routine called by sbcmod_tam
75   PUBLIC   sbc_fwb_adj_tst  ! routine called by tst
76
77   REAL(wp), PUBLIC ::   a_fwb_tl           ! for before year.
78   REAL(wp), PUBLIC ::   a_fwb_ad           ! for before year.
79   REAL(wp) ::   area               ! global mean ocean surface (interior domain)
80
81   REAL(wp), DIMENSION(jpi,jpj) ::  &
82     & e1e2_i    ! area of the interior domain (e1t*e2t*tmask_i)
83
84   !! * Substitutions
85#  include "domzgr_substitute.h90"
86#  include "vectopt_loop_substitute.h90"
87   !!----------------------------------------------------------------------
88   !!  OPA 9.0 , LOCEAN-IPSL (2006)
89   !! $Id: sbcfwb.F90 1168 2008-08-11 10:21:06Z rblod $
90   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
91   !!----------------------------------------------------------------------
92CONTAINS
93
94   SUBROUTINE sbc_fwb_tan( kt, kn_fwb, kn_fsbc )
95      !!---------------------------------------------------------------------
96      !!                  ***  ROUTINE sbc_fwb_tan  ***
97      !!
98      !! ** Purpose :   Control the mean sea surface drift
99      !!
100      !! ** Method  :   several ways  depending on kn_fwb
101      !!                =0 no control
102      !!                =1 global mean of emp set to zero at each nn_fsbc time step
103      !!                =2 annual global mean corrected from previous year
104      !!                =3 variable relaxation time to zero or to a time varying field
105      !!----------------------------------------------------------------------
106      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
107      INTEGER, INTENT( in ) ::   kn_fsbc  !
108      INTEGER, INTENT( in ) ::   kn_fwb   ! ocean time-step index
109      !!
110      INTEGER  ::   inum                  ! temporary logical unit
111      INTEGER  ::   ikty, iyear           !
112      CHARACTER (len=32) ::   clname 
113      REAL(wp) ::   z_emptl               ! temporary scalars
114      !!----------------------------------------------------------------------
115      !
116      IF( kt == nit000 ) THEN
117         !
118         IF(lwp) THEN
119            WRITE(numout,*)
120            WRITE(numout,*) 'sbc_fwb_tan : FreshWater Budget correction.'
121            WRITE(numout,*) '~~~~~~~~~~~'
122         ENDIF
123         !
124         e1e2_i(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
125         area = SUM( e1e2_i(:,:) )
126         IF( lk_mpp )   CALL  mpp_sum( area    )   ! sum over the global domain
127         !
128      ENDIF
129     
130
131      SELECT CASE ( kn_fwb )
132      !
133      CASE ( 1 )                               ! global mean emp set to zero
134         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
135            z_emptl = SUM( e1e2_i(:,:) * emp_tl(:,:) ) / area
136            IF( lk_mpp )   CALL  mpp_sum( z_emptl    )   ! sum over the global domain
137            emp_tl (:,:) = emp_tl (:,:) - z_emptl
138            emps_tl(:,:) = emps_tl(:,:) - z_emptl
139         ENDIF
140         !
141      CASE ( 2 )                               ! emp budget adjusted from the previous year
142         ! initialisation
143         IF( kt == nit000 ) THEN
144            a_fwb_tl = 0.0_wp
145         ENDIF   
146         !
147         ikty = 365 * 86400 / rdttra(1)    !!bug  use of 365 days leap year or 360d year !!!!!!!
148         IF( MOD( kt, ikty ) == 0 ) THEN
149            a_fwb_tl   = SUM( e1e2_i(:,:) * sshn_tl(:,:) )
150            IF( lk_mpp )   CALL  mpp_sum( a_fwb_tl )   ! sum over the global domain
151            a_fwb_tl   = a_fwb_tl * 1.e+3 / ( area * 86400. * 365. )     ! convert in Kg/m3/s = mm/s
152!!gm        !                                                      !!bug 365d year
153         ENDIF
154         !
155         ! correct the freshwater fluxes
156         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
157            emp_tl (:,:) = emp_tl (:,:) + a_fwb_tl
158            emps_tl(:,:) = emps_tl(:,:) + a_fwb_tl
159         ENDIF
160         !
161      CASE ( 3 )
162         !
163         WRITE(ctmp1,*)'sbc_fwb_tan: nn_fwb=', kn_fwb, ' is not available for TAM yet , choose either 0/1/2'
164         CALL ctl_stop( ctmp1 )
165
166         !
167      CASE DEFAULT    ! you should never be there
168         WRITE(ctmp1,*)'sbc_fwb_tan: nn_fwb=', kn_fwb, ' is not permitted for the FreshWater Budget correction, choose either 0/1/2'
169         CALL ctl_stop( ctmp1 )
170         !
171      END SELECT
172      !
173   END SUBROUTINE sbc_fwb_tan
174
175
176   SUBROUTINE sbc_fwb_adj( kt, kn_fwb, kn_fsbc )
177      !!---------------------------------------------------------------------
178      !!                  ***  ROUTINE sbc_fwb_adj  ***
179      !!
180      !! ** Purpose :   Control the mean sea surface drift
181      !!
182      !! ** Method  :   several ways  depending on kn_fwb
183      !!                =0 no control
184      !!                =1 global mean of emp set to zero at each nn_fsbc time step
185      !!                =2 annual global mean corrected from previous year
186      !!                =3 variable relaxation time to zero or to a time varying field
187      !!----------------------------------------------------------------------
188      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
189      INTEGER, INTENT( in ) ::   kn_fsbc  !
190      INTEGER, INTENT( in ) ::   kn_fwb   ! ocean time-step index
191      !!
192      INTEGER  ::   inum                  ! temporary logical unit
193      INTEGER  ::   ikty, iyear           !
194      INTEGER  ::   ji, jj           !
195      CHARACTER (len=32) ::   clname 
196      REAL(wp) ::   z_empad               ! temporary scalars
197      !!----------------------------------------------------------------------
198      !
199      IF( kt == nitend ) THEN
200         !
201         IF(lwp) THEN
202            WRITE(numout,*)
203            WRITE(numout,*) 'sbc_fwb_adj : FreshWater Budget correction.'
204            WRITE(numout,*) '~~~~~~~~~~~'
205         ENDIF
206         !
207         e1e2_i(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
208         area = SUM( e1e2_i(:,:) )
209         IF( lk_mpp )   CALL  mpp_sum( area    )   ! sum over the global domain
210         !
211         ! initialisation
212         a_fwb_ad = 0.0_wp
213      ENDIF
214         z_empad  = 0.0_wp     
215
216      SELECT CASE ( kn_fwb )
217      CASE ( 0 )
218         ! No Control
219      !
220      CASE ( 1 )                               ! global mean emp set to zero
221         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
222            DO jj = 1, jpj
223               DO ji = 1, jpi
224                  z_empad = z_empad - emps_ad(ji,jj) - emp_ad(ji,jj)
225               END DO
226            END DO
227            IF( lk_mpp )   CALL  mpp_sum( z_empad    )   ! sum over the global domain
228            z_empad = z_empad / area
229            DO jj = 1, jpj
230               DO ji = 1, jpi
231                  emp_ad(ji,jj) = emp_ad(ji,jj) +  e1e2_i(ji,jj) * z_empad
232               END DO
233            END DO
234         ENDIF
235         !
236      CASE ( 2 )                               ! emp budget adjusted from the previous year
237         !
238         ! initialisation
239         ikty = 365 * 86400 / rdttra(1)    !!bug  use of 365 days leap year or 360d year !!!!!!!
240         ! correct the freshwater fluxes
241         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
242            DO jj = nldj, nlej
243               DO ji = nldi, nlei
244                  a_fwb_ad = a_fwb_ad + emps_ad(ji,jj) + emp_ad (ji,jj)
245               END DO
246            END DO
247         ENDIF
248         !
249         IF( MOD( kt, ikty ) == 0 ) THEN
250            a_fwb_ad   = a_fwb_ad * 1.e+3 / ( area * 86400. * 365. )     ! convert in Kg/m3/s = mm/s
251            IF( lk_mpp )   CALL  mpp_sum( a_fwb_ad )   ! sum over the global domain
252            DO jj = 1, jpj
253               DO ji = 1, jpi
254                  sshn_ad(ji,jj) = sshn_ad(ji,jj) + e1e2_i(ji,jj) * a_fwb_ad
255               END DO
256            END DO
257            a_fwb_ad = 0.0_wp
258         ENDIF   
259         !
260      CASE ( 3 )
261         !
262         WRITE(ctmp1,*)'sbc_fwb_tan: nn_fwb=', kn_fwb, ' is not available for TAM yet , choose either 0/1/2'
263         CALL ctl_stop( ctmp1 )
264         !
265      CASE DEFAULT    ! you should never be there
266         WRITE(ctmp1,*)'sbc_fwb_adj: nn_fwb=', kn_fwb, ' is not permitted for the FreshWater Budget correction, choose either 0/1/2'
267         CALL ctl_stop( ctmp1 )
268         !
269      END SELECT
270      !
271   END SUBROUTINE sbc_fwb_adj
272
273   SUBROUTINE sbc_fwb_adj_tst ( kumadt )
274      !!-----------------------------------------------------------------------
275      !!
276      !!                  ***  ROUTINE example_adj_tst ***
277      !!
278      !! ** Purpose : Test the adjoint routine.
279      !!
280      !! ** Method  : Verify the scalar product
281      !!           
282      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
283      !!
284      !!              where  L   = tangent routine
285      !!                     L^T = adjoint routine
286      !!                     W   = diagonal matrix of scale factors
287      !!                     dx  = input perturbation (random field)
288      !!                     dy  = L dx
289      !!                   
290      !! History :
291      !!        ! 08-08 (A. Vidard)
292      !!-----------------------------------------------------------------------
293      !! * Modules used
294
295      !! * Arguments
296      INTEGER, INTENT(IN) :: &
297         & kumadt             ! Output unit
298 
299      !! * Local declarations
300      INTEGER ::  &
301         & istp,  &
302         & jstp,  &
303         & ji,    &        ! dummy loop indices
304         & jj,    &       
305         & jk,    &
306         & jn_fwb   
307      INTEGER, DIMENSION(jpi,jpj) :: &
308         & iseed_2d        ! 2D seed for the random number generator
309      REAL(KIND=wp) :: &
310         & zsp1,         & ! scalar product involving the tangent routine
311         & zsp2            ! scalar product involving the adjoint routine
312      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
313         & zemp_tlin ,     & ! Tangent input
314         & zemps_tlin ,    & ! Tangent input
315         & zssh_tlin ,     & ! Tangent input
316#if defined key_obc
317         & zssh_tlout, zssh_adin, &
318#endif
319         & zemp_tlout,     & ! Tangent output
320         & zemps_tlout,    & ! Tangent output
321         & zemp_adin ,     & ! Adjoint input
322         & zemps_adin ,    & ! Adjoint input
323         & zemp_adout,     & ! Adjoint output
324         & zemps_adout,    & ! Adjoint output
325         & zssh_adout,     & ! Adjoint output
326         & zr                ! 3D random field
327      CHARACTER(LEN=14) :: cl_name
328      ! Allocate memory
329
330      ALLOCATE( &
331#if defined key_obc
332         & zssh_tlout( jpi,jpj),  zssh_adin ( jpi,jpj), &
333#endif
334         & zemp_tlin(  jpi,jpj),     &
335         & zemps_tlin( jpi,jpj),     &
336         & zssh_tlin(  jpi,jpj),     &
337         & zemp_tlout( jpi,jpj),     &
338         & zemps_tlout(jpi,jpj),     &
339         & zemp_adin(  jpi,jpj),     &
340         & zemps_adin( jpi,jpj),     &
341         & zemp_adout( jpi,jpj),     &
342         & zemps_adout(jpi,jpj),     &
343         & zssh_adout( jpi,jpj),     &
344         & zr(         jpi,jpj)      &
345         & )
346      !==================================================================
347      ! 1) dx = ( emp_tl, emps_tl, ssh_tl ) and
348      !    dy = ( emp_tl, emps_tl )
349      !==================================================================
350      ! Test for time steps nit000 and nit000 + 1 (the matrix changes)
351
352      DO jstp = nit000, nit000 + 2
353         istp = jstp
354         IF ( jstp == nit000 +2 ) istp = nitend
355
356      !--------------------------------------------------------------------
357      ! Reset the tangent and adjoint variables
358      !--------------------------------------------------------------------
359          zemp_tlin  (:,:) = 0.0_wp
360          zemps_tlin (:,:) = 0.0_wp
361          zssh_tlin  (:,:) = 0.0_wp
362          zemp_tlout (:,:) = 0.0_wp
363          zemps_tlout(:,:) = 0.0_wp
364          zemp_adin  (:,:) = 0.0_wp
365          zemps_adin (:,:) = 0.0_wp
366          zemp_adout (:,:) = 0.0_wp
367          zemps_adout(:,:) = 0.0_wp
368          zssh_adout (:,:) = 0.0_wp
369          zr(:,:)          = 0.0_wp
370
371      !--------------------------------------------------------------------
372      ! Initialize the tangent input with random noise: dx
373      !--------------------------------------------------------------------
374
375      DO jj = 1, jpj
376         DO ji = 1, jpi
377            iseed_2d(ji,jj) = - ( 596035 + &
378               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
379         END DO
380      END DO
381      CALL grid_random( iseed_2d, zr, 'T', 0.0_wp, stdemp )
382      DO jj = nldj, nlej
383         DO ji = nldi, nlei
384            zemp_tlin(ji,jj) = zr(ji,jj) 
385         END DO
386      END DO
387      DO jj = 1, jpj
388         DO ji = 1, jpi
389            iseed_2d(ji,jj) = - ( 446251 + &
390               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
391         END DO
392      END DO
393      CALL grid_random( iseed_2d, zr, 'T', 0.0_wp, stdemp )
394      DO jj = nldj, nlej
395         DO ji = nldi, nlei
396            zemps_tlin(ji,jj) = zr(ji,jj) 
397         END DO
398      END DO 
399
400      DO jj = 1, jpj
401         DO ji = 1, jpi
402            iseed_2d(ji,jj) = - ( 352678 + &
403               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
404         END DO
405      END DO
406      CALL grid_random( iseed_2d, zr, 'T', 0.0_wp, stdssh )
407      DO jj = nldj, nlej
408         DO ji = nldi, nlei
409            zssh_tlin(ji,jj) = zr(ji,jj) 
410         END DO
411      END DO
412
413      DO jn_fwb = 1, 2
414
415         a_fwb_tl         = 0.0_wp
416         a_fwb_ad         = 0.0_wp
417         sshn_ad    (:,:) = 0.0_wp
418
419         sshn_tl(:,:) = zssh_tlin (:,:)
420         emps_tl(:,:) = zemps_tlin(:,:)
421         emp_tl (:,:) = zemp_tlin (:,:)
422
423         CALL sbc_fwb_tan( istp, jn_fwb, 1 )
424
425         zemps_tlout(:,:) = emps_tl(:,:)
426         zemp_tlout (:,:) = emp_tl (:,:)
427#if defined key_obc
428         zssh_tlout (:,:) = sshn_tl (:,:)
429#endif
430         !-----------------------------------------------------------------
431         ! Initialize the adjoint variables: dy^* = W dy
432         !-----------------------------------------------------------------
433
434         DO jj = nldj, nlej
435            DO ji = nldi, nlei
436               zemp_adin( ji,jj) = zemp_tlout( ji,jj) &
437                  &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) &
438                  &               * tmask(ji,jj,1)
439               zemps_adin(ji,jj) = zemps_tlout(ji,jj) &
440                  &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) &
441                  &               * tmask(ji,jj,1)
442#if defined key_obc
443               zssh_adin( ji,jj) = zssh_tlout( ji,jj) * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) * tmask(ji,jj,1)
444#endif
445            END DO
446         END DO
447
448         !-----------------------------------------------------------------
449         ! Compute the scalar product: ( L dx )^T W dy
450         !-----------------------------------------------------------------
451
452         zsp1 = DOT_PRODUCT( zemp_tlout,  zemp_adin  ) &
453            & + DOT_PRODUCT( zemps_tlout, zemps_adin )
454
455#if defined key_obc
456         zsp1 = zsp1 + DOT_PRODUCT( zssh_tlout, zssh_adin )
457#endif
458         !-----------------------------------------------------------------
459         ! Call the adjoint routine: dx^* = L^T dy^*
460         !-----------------------------------------------------------------
461
462         emp_ad (:,:) = zemp_adin (:,:)
463         emps_ad(:,:) = zemps_adin(:,:)
464#if defined key_obc
465         sshn_ad(:,:) = zssh_adin(:,:)
466#endif
467         CALL sbc_fwb_adj ( istp, jn_fwb, 1 )
468
469         zemps_adout(:,:) = emps_ad(:,:)
470         zemp_adout (:,:) = emp_ad (:,:)
471         zssh_adout (:,:) = sshn_ad(:,:)
472
473         zsp2 = DOT_PRODUCT( zemp_tlin,  zemp_adout  ) &
474            & + DOT_PRODUCT( zssh_tlin,  zssh_adout  ) &
475            & + DOT_PRODUCT( zemps_tlin, zemps_adout ) 
476
477         ! 14 char:'12345678901234'
478         IF ( istp == nit000 ) THEN
479             WRITE (cl_name,"(A10,1x,i1,1x,A1)") 'sbcfwb_adj',jn_fwb,'1'
480         ELSEIF ( istp == nit000 + 1 ) THEN
481             WRITE (cl_name,"(A10,1x,i1,1x,A1)") 'sbcfwb_adj',jn_fwb,'2'
482         ELSEIF ( istp == nitend ) THEN
483             WRITE (cl_name,"(A10,1x,i1,1x,A1)") 'sbcfwb_adj',jn_fwb,'3'
484         END IF
485!         WRITE (cl_name,"(A11,2x,i1)") 'sbc_fwb_adj',jn_fwb
486         CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
487
488      END DO
489
490      END DO
491
492      DEALLOCATE(       &
493         & zemp_tlin,   &
494         & zemps_tlin,  &
495         & zssh_tlin,   &
496         & zemp_tlout,  &
497         & zemps_tlout, &
498         & zemp_adin,   &
499         & zemps_adin,  &
500         & zemp_adout,  &
501         & zemps_adout, &
502         & zssh_adout,  &
503         & zr           &
504         & )
505
506   END SUBROUTINE sbc_fwb_adj_tst
507#endif
508   !!======================================================================
509END MODULE sbcfwb_tam
Note: See TracBrowser for help on using the repository browser.