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/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/SBC – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/SBC/sbcfwb_tam.F90 @ 3611

Last change on this file since 3611 was 3611, checked in by pabouttier, 11 years ago

Add TAM code and ORCA2_TAM configuration - see Ticket #1007

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