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.
tstool_tam.F90 in NEMO/branches/NERC/dev_release-3.4_NEMOTAM_consolidated/NEMOGCM/NEMO/OPATAM_SRC – NEMO

source: NEMO/branches/NERC/dev_release-3.4_NEMOTAM_consolidated/NEMOGCM/NEMO/OPATAM_SRC/tstool_tam.F90 @ 13432

Last change on this file since 13432 was 5168, checked in by pabouttier, 9 years ago

Fixed several bugs described in Ticket #1360

  • Property svn:executable set to *
File size: 7.8 KB
Line 
1MODULE tstool_tam
2   !!==========================================================================
3   !!     ***  MODULE  tstool_tam : TAM testing utilities  ***
4   !!==========================================================================
5   !! History of the NEMOTAM module:
6   !!             3.0  !  08-11  (A. Vidard) initial version
7   USE par_oce,        ONLY: & ! Ocean space and time domain variables
8      & jpi,                 &
9      & jpj,                 &
10      & jpk,                 &
11      & jpiglo,              &
12      & jpjglo
13   USE dom_oce       , ONLY: & ! Ocean space and time domain
14      & e1u,                 &
15      & e2u,                 &
16      & e3t_0,               &
17      & e3u,                 &
18      & umask,               &
19      & mig,                 &
20      & mjg,                 &
21      & nldi,                &
22      & nldj,                &
23      & nlei,                &
24      & nlej
25   USE wrk_nemo       ! Memory Allocation
26   USE par_kind      , ONLY: & ! Precision variables
27      & wp
28   USE in_out_manager, ONLY: & ! I/O manager
29      & lwp
30   USE gridrandom    , ONLY: & ! Random Gaussian noise on grids
31      & grid_random
32   USE dotprodfld    , ONLY: & ! Computes dot product for 3D and 2D fields
33      & dot_product
34   IMPLICIT NONE
35   PRIVATE
36   REAL(KIND=wp), PUBLIC ::  & ! random field standard deviation for:
37      & stdu   =  0.01_wp,    & !   u-velocity
38      & stdv   =  0.01_wp,    & !   v-velocity
39      & stdw   = 0.01_wp,    & !   w-velocity
40#if defined key_obc
41      & stds   = 0.001_wp,    & !   salinity
42      & stdt   = 0.020_wp,    & !   temperature
43      & stdssh = 0.005_wp,   & !   sea surface height
44#else
45      & stds   =  0.01_wp,    & !   salinity
46      & stdt   =  0.1_wp,    & !   temperature
47      & stdssh = 0.01_wp,    & !   sea surface height
48#endif
49      & stdemp = 0.001_wp,    & !   evaporation minus precip 0.1_wp / SQRT( wesp_emp )
50      & stdqns =  0.1_wp,    & !   non solar heat flux
51      & stdqsr =  0.1_wp,    & !   solar heat flux
52      & stdgc  =  0.01_wp,    & !   gcx, gcb
53      & stdr   =  0.01_wp,    & !   rotb, rhd
54      & stdh   =  0.01_wp       !   hdivb
55
56   PUBLIC &
57      & prntst_adj,          &
58      & prntst_tlm
59
60#  include "domzgr_substitute.h90"
61
62CONTAINS
63
64   SUBROUTINE prntst_adj( cd_name, kumadt, psp1, psp2 )
65      CHARACTER(LEN=14), INTENT(in) :: cd_name
66      REAL(wp), INTENT(in) :: psp1, psp2
67      INTEGER, INTENT(in) :: kumadt
68      REAL(KIND=wp) :: &
69         & zspdif,       & ! scalar product difference
70         & ztol            ! accepted tolerance
71      CHARACTER (LEN=47) :: &
72         & FMT
73      CHARACTER (LEN=9)  :: &
74         & cl_stat         ! Accuracy status of adjoint routine (ok or warning)
75
76      ! Compare the scalar products
77
78      zspdif = ABS( psp1 - psp2 )
79      IF ( psp1 /= 0.0_wp ) zspdif = zspdif / ABS( psp1 )
80
81      ztol = EPSILON( zspdif ) * 10.
82
83      IF ( zspdif < ztol ) THEN
84         cl_stat = '   ok    '
85      ELSEIF ( zspdif < ztol*1000._wp ) THEN
86         cl_stat = ' warning '
87      ELSE
88         cl_stat = 'RED ALERT'
89      ENDIF
90
91      IF (lwp) THEN
92         FMT = "(A14,1X,E20.15,2X,E20.15,2X,E6.1,1X,E6.1,1x,A9)"
93         WRITE(kumadt,FMT) cd_name, psp1, psp2, zspdif, ztol, cl_stat
94         CALL FLUSH( kumadt )
95      ENDIF
96   END SUBROUTINE prntst_adj
97
98   SUBROUTINE prntst_tlm( cd_name, kumadt, psp1, psp2 )
99      CHARACTER(LEN=14), INTENT(in) :: cd_name
100      REAL(wp), INTENT(in) :: psp1, psp2
101      INTEGER, INTENT(in) :: kumadt
102      REAL(KIND=wp) :: &
103         & zspratio       ! scalar product difference
104      CHARACTER (LEN=47) :: &
105         & FMT
106      CHARACTER (LEN=9)  :: &
107         & cl_stat         ! Accuracy status of adjoint routine (ok or warning)
108
109      ! Compare the scalar products
110
111      IF ( psp1 /= 0.0_wp ) zspratio = 100 * psp1 /  psp2
112
113      IF (lwp) THEN
114         FMT = "(A14,1X,E20.13,2X,E20.15,2X,E6.1,1X)"
115         WRITE(kumadt,FMT) cd_name, psp1, psp2, zspratio
116         CALL FLUSH( kumadt )
117      ENDIF
118   END SUBROUTINE prntst_tlm
119
120   SUBROUTINE example_adj_tst( kumadt )
121      !!-----------------------------------------------------------------------
122      !!
123      !!                  ***  ROUTINE example_adj_tst ***
124      !!
125      !! ** Purpose : Test the adjoint routine.
126      !!
127      !! ** Method  : Verify the scalar product
128      !!
129      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
130      !!
131      !!              where  L   = tangent routine
132      !!                     L^T = adjoint routine
133      !!                     W   = diagonal matrix of scale factors
134      !!                     dx  = input perturbation (random field)
135      !!                     dy  = L dx
136      !!
137      !!
138      !! History :
139      !!        ! 08-08 (A. Vidard)
140      !!-----------------------------------------------------------------------
141      !! * Modules used
142
143      !! * Arguments
144      INTEGER, INTENT(IN) :: &
145         & kumadt             ! Output unit
146
147      !! * Local declarations
148      INTEGER ::  &
149         & ji,    &        ! dummy loop indices
150         & jj,    &
151         & jk
152      REAL(KIND=wp) :: &
153         & zsp1,         & ! scalar product involving the tangent routine
154         & zsp2            ! scalar product involving the adjoint routine
155      REAL(KIND=wp), POINTER, DIMENSION(:,:,:) :: &
156         & z_tlin ,     & ! Tangent input
157         & z_tlout,     & ! Tangent output
158         & z_adin ,     & ! Adjoint input
159         & z_adout,     & ! Adjoint output
160         & zr             ! 3D random field
161      CHARACTER(LEN=14) :: &
162         & cl_name
163      ! Allocate memory
164
165      CALL wrk_alloc( jpi, jpj, jpk, z_tlin , z_tlout, z_adin, z_adout, zr )
166
167      !==================================================================
168      ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
169      !    dy = ( hdivb_tl, hdivn_tl )
170      !==================================================================
171
172      !--------------------------------------------------------------------
173      ! Reset the tangent and adjoint variables
174      !--------------------------------------------------------------------
175          z_tlin( :,:,:) = 0.0_wp
176          z_tlout(:,:,:) = 0.0_wp
177          z_adin( :,:,:) = 0.0_wp
178          z_adout(:,:,:) = 0.0_wp
179          zr(     :,:,:) = 0.0_wp
180      !--------------------------------------------------------------------
181      ! Initialize the tangent input with random noise: dx
182      !--------------------------------------------------------------------
183
184      CALL grid_random( zr, 'U', 0.0_wp, stdr )
185      z_tlin(:,:,:) = zr(:,:,:)
186
187      CALL example_tan
188      !--------------------------------------------------------------------
189      ! Initialize the adjoint variables: dy^* = W dy
190      !--------------------------------------------------------------------
191
192      DO jk = 1, jpk
193        DO jj = nldj, nlej
194           DO ji = nldi, nlei
195              z_adin(ji,jj,jk) = z_tlout(ji,jj,jk) &
196                 &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
197                 &               * umask(ji,jj,jk)
198            END DO
199         END DO
200      END DO
201      !--------------------------------------------------------------------
202      ! Compute the scalar product: ( L dx )^T W dy
203      !--------------------------------------------------------------------
204
205      zsp1 = DOT_PRODUCT( z_tlout, z_adin )
206
207      !--------------------------------------------------------------------
208      ! Call the adjoint routine: dx^* = L^T dy^*
209      !--------------------------------------------------------------------
210
211      CALL example_adj
212
213      zsp2 = DOT_PRODUCT( z_tlin, z_adout )
214
215      ! 14 char:'12345678901234'
216      cl_name = 'example_adj   '
217      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
218
219      CALL wrk_dealloc( jpi, jpj, jpk, z_tlin , z_tlout, z_adin, z_adout, zr )
220
221   END SUBROUTINE example_adj_tst
222
223   SUBROUTINE example_tan
224   END SUBROUTINE example_tan
225   SUBROUTINE example_adj
226   END SUBROUTINE example_adj
227
228END MODULE tstool_tam
Note: See TracBrowser for help on using the repository browser.