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.
closea_tam.F90 in branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DOM – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DOM/closea_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: 18.2 KB
Line 
1MODULE closea_tam
2#if defined key_tam
3   !!======================================================================
4   !!                       ***  MODULE  closea_tam  ***
5   !! Closed Seas  : specific treatments associated with closed seas
6   !!                Tangent and adjoint module
7   !!======================================================================
8   !! History of the direct module:
9   !!             8.2  !  00-05  (O. Marti)  Original code
10   !!             8.5  !  02-06  (E. Durand, G. Madec)  F90
11   !!             9.0  !  06-07  (G. Madec)  add clo_rnf, clo_ups, clo_bat
12   !! History of the tangent module:
13   !!             9.0  !  08-11  (A. Vidard)   skeleton tam of sbc_clo
14   !!             9.0  !  09-08  (F. Vigilant) tangent and adjoint version
15   !!        NEMO 3.4  !  07-12  (P.-A. Bouttier) Phasing with 3.4
16   !!----------------------------------------------------------------------
17
18   !!----------------------------------------------------------------------
19   !!   sbc_clo    : Special handling of closed seas
20   !!----------------------------------------------------------------------
21
22   USE par_oce
23   USE par_kind
24   USE sbc_oce_tam
25   USE lbclnk_tam
26   USE lbclnk
27   USE dom_oce
28   USE sbc_oce
29   USE closea
30   USE lib_mpp                 ! distribued memory computing library
31   USE in_out_manager
32   USE gridrandom
33   USE dotprodfld
34   USE paresp
35   USE tstool_tam
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC sbc_clo_tan      ! routine called by sbcmod_tam module
41   PUBLIC sbc_clo_adj      ! routine called by sbcmod_tam module
42   PUBLIC sbc_clo_adj_tst  ! routine called by tst module
43
44   REAL(wp), DIMENSION (jpncs+1)       ::   surf             ! closed sea surface
45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !!  OPA 9.0 , LOCEAN-IPSL (2006)
51   !! $Id: closea.F90 1146 2008-06-25 11:42:56Z rblod $
52   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54
55CONTAINS
56   SUBROUTINE sbc_clo_tan( kt )
57      !!---------------------------------------------------------------------
58      !!                  ***  ROUTINE sbc_clo_tan  ***
59      !!
60      !! ** Purpose :   Special handling of closed seas (Tangent version)
61      !!
62      !! ** Method  :   Water flux is forced to zero over closed sea
63      !!      Excess is shared between remaining ocean, or
64      !!      put as run-off in open ocean.
65      !!
66      !! ** Action  :   emp, emps   updated surface freshwater fluxes at kt
67      !!----------------------------------------------------------------------
68      INTEGER, INTENT(in) ::   kt   ! ocean model time step
69      !
70      INTEGER                     ::   ji, jj, jc, jn   ! dummy loop indices
71      REAL(wp)                    ::   zze2
72      REAL(wp), DIMENSION (jpncs) ::   zfwftl
73      !!----------------------------------------------------------------------
74      !
75      !                                                   !------------------!
76      IF( kt == nit000 ) THEN                             !  Initialisation  !
77         !                                                !------------------!
78         IF(lwp) WRITE(numout,*)
79         IF(lwp) WRITE(numout,*)'sbc_clo_tan : closed seas '
80         IF(lwp) WRITE(numout,*)'~~~~~~~'
81
82         ! Total surface of ocean
83         surf(jpncs+1) = SUM( e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
84
85         DO jc = 1, jpncs
86            surf(jc) =0.0_wp
87            DO jj = ncsj1(jc), ncsj2(jc)
88               DO ji = ncsi1(jc), ncsi2(jc)
89                  surf(jc) = surf(jc) + e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)      ! surface of closed seas
90               END DO
91            END DO
92         END DO
93         IF( lk_mpp )   CALL mpp_sum ( surf, jpncs+1 )       ! mpp: sum over all the global domain
94
95         IF(lwp) WRITE(numout,*)'     Closed sea surfaces'
96         DO jc = 1, jpncs
97            IF(lwp)WRITE(numout,FMT='(1I3,4I4,5X,F16.2)') jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc)
98         END DO
99
100         ! jpncs+1 : surface of sea, closed seas excluded
101         DO jc = 1, jpncs
102            surf(jpncs+1) = surf(jpncs+1) - surf(jc)
103         END DO
104         !
105      ENDIF
106      !                                                   !--------------------!
107      !                                                   !  update emp, emps  !
108      zfwftl = 0.0_wp                                     !--------------------!
109      DO jc = 1, jpncs
110         DO jj = ncsj1(jc), ncsj2(jc)
111            DO ji = ncsi1(jc), ncsi2(jc)
112               zfwftl(jc) = zfwftl(jc) + e1t(ji,jj) * e2t(ji,jj) * ( emp_tl(ji,jj) - rnf(ji,jj) ) * tmask_i(ji,jj)
113            END DO
114         END DO
115      END DO
116      IF( lk_mpp )   CALL mpp_sum ( zfwftl(:) , jpncs )       ! mpp: sum over all the global domain
117
118      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN      ! Black Sea case for ORCA_R2 configuration
119         zze2    = ( zfwftl(3) + zfwftl(4) ) / 2.
120         zfwftl(3) = zze2
121         zfwftl(4) = zze2
122      ENDIF
123      DO jc = 1, jpncs
124         !
125         IF( ncstt(jc) == 0 ) THEN
126            ! water/evap excess is shared by all open ocean
127            emp_tl (:,:) = emp_tl (:,:) + zfwftl(jc) / surf(jpncs+1)
128            emps_tl(:,:) = emps_tl(:,:) + zfwftl(jc) / surf(jpncs+1)
129         ELSEIF( ncstt(jc) == 1 ) THEN
130            ! Excess water in open sea, at outflow location, excess evap shared
131            IF ( zfwftl(jc) <= 0.e0 ) THEN
132                DO jn = 1, ncsnr(jc)
133                  ji = mi0(ncsir(jc,jn))
134                  jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
135                  IF (      ji > 1 .AND. ji < jpi   &
136                      .AND. jj > 1 .AND. jj < jpj ) THEN
137                      emp_tl (ji,jj) = emp_tl (ji,jj) + zfwftl(jc) /   &
138                         (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj))
139                      emps_tl(ji,jj) = emps_tl(ji,jj) + zfwftl(jc) /   &
140                          (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj))
141                  END IF
142                END DO
143            ELSE
144                emp_tl (:,:) = emp_tl (:,:) + zfwftl(jc) / surf(jpncs+1)
145                emps_tl(:,:) = emps_tl(:,:) + zfwftl(jc) / surf(jpncs+1)
146            ENDIF
147         ELSEIF( ncstt(jc) == 2 ) THEN
148            ! Excess e-p+r (either sign) goes to open ocean, at outflow location
149            IF(      ji > 1 .AND. ji < jpi    &
150               .AND. jj > 1 .AND. jj < jpj ) THEN
151                DO jn = 1, ncsnr(jc)
152                  ji = mi0(ncsir(jc,jn))
153                  jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
154                  emp_tl (ji,jj) = emp_tl (ji,jj) + zfwftl(jc)   &
155                      / (FLOAT(ncsnr(jc)) *  e1t(ji,jj) * e2t(ji,jj) )
156                  emps_tl(ji,jj) = emps_tl(ji,jj) + zfwftl(jc)   &
157                      / (FLOAT(ncsnr(jc)) *  e1t(ji,jj) * e2t(ji,jj) )
158                END DO
159            ENDIF
160         ENDIF
161         !
162         DO jj = ncsj1(jc), ncsj2(jc)
163            DO ji = ncsi1(jc), ncsi2(jc)
164               emp_tl (ji,jj) = emp_tl (ji,jj) - zfwftl(jc) / surf(jc)
165               emps_tl(ji,jj) = emps_tl(ji,jj) - zfwftl(jc) / surf(jc)
166            END DO
167         END DO
168         !
169      END DO
170      !
171      CALL lbc_lnk( emp_tl , 'T', 1. )
172      CALL lbc_lnk( emps_tl, 'T', 1. )
173      !
174   END SUBROUTINE sbc_clo_tan
175   SUBROUTINE sbc_clo_adj( kt )
176      !!---------------------------------------------------------------------
177      !!                  ***  ROUTINE sbc_clo_adj  ***
178      !!
179      !! ** Purpose :   Special handling of closed seas (Adjoint version)
180      !!
181      !! ** Method  :   Water flux is forced to zero over closed sea
182      !!      Excess is shared between remaining ocean, or
183      !!      put as run-off in open ocean.
184      !!
185      !! ** Action  :   emp, emps   updated surface freshwater fluxes at kt
186      !!----------------------------------------------------------------------
187      INTEGER, INTENT(in) ::   kt   ! ocean model time step
188      !
189      !
190      INTEGER                     ::   ji, jj, jc, jn   ! dummy loop indices
191      REAL(wp)                    ::   zze2
192      REAL(wp), DIMENSION (jpncs) ::   zfwfad
193      !!----------------------------------------------------------------------
194      !                                                   !------------------!
195      IF( kt == nit000 ) THEN                             !  Initialisation  !
196         !                                                !------------------!
197         IF(lwp) WRITE(numout,*)
198         IF(lwp) WRITE(numout,*)'sbc_clo_adj : closed seas '
199         IF(lwp) WRITE(numout,*)'~~~~~~~'
200
201         ! Total surface of ocean
202         surf(jpncs+1) = SUM( e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
203
204         DO jc = 1, jpncs
205            surf(jc) =0.0_wp
206            DO jj = ncsj1(jc), ncsj2(jc)
207               DO ji = ncsi1(jc), ncsi2(jc)
208                  surf(jc) = surf(jc) + e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)      ! surface of closed seas
209               END DO
210            END DO
211         END DO
212         IF( lk_mpp )   CALL mpp_sum ( surf, jpncs+1 )       ! mpp: sum over all the global domain
213
214         IF(lwp) WRITE(numout,*)'     Closed sea surfaces'
215         DO jc = 1, jpncs
216            IF(lwp)WRITE(numout,FMT='(1I3,4I4,5X,F16.2)') jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc)
217         END DO
218
219         ! jpncs+1 : surface of sea, closed seas excluded
220         DO jc = 1, jpncs
221            surf(jpncs+1) = surf(jpncs+1) - surf(jc)
222         END DO
223         !
224      ENDIF
225      zfwfad = 0.0_wp
226      CALL lbc_lnk_adj( emp_ad , 'T', 1. )
227      CALL lbc_lnk_adj( emps_ad, 'T', 1. )
228      DO jc = jpncs, 1, -1
229         !
230         DO jj = ncsj1(jc), ncsj2(jc)
231            DO ji = ncsi1(jc), ncsi2(jc)
232               zfwfad(jc) = zfwfad(jc) - emps_ad(ji,jj) / surf(jc)
233               zfwfad(jc) = zfwfad(jc) - emp_ad(ji,jj)  / surf(jc)
234            END DO
235         END DO
236         !
237         IF( ncstt(jc) == 0 ) THEN
238            ! water/evap excess is shared by all open ocean
239            DO jj = 1, jpj
240               DO ji = 1, jpi
241                  zfwfad(jc) = zfwfad(jc) + emps_ad(ji,jj) / surf(jpncs+1)
242                  zfwfad(jc) = zfwfad(jc) + emp_ad(ji,jj)  / surf(jpncs+1)
243               END DO
244            END DO
245         ELSEIF( ncstt(jc) == 1 ) THEN
246            ! Excess water in open sea, at outflow location, excess evap shared
247            IF ( zfwfad(jc) <= 0.e0 ) THEN
248                DO jn = 1, ncsnr(jc)
249                  ji = mi0(ncsir(jc,jn))
250                  jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
251                  IF (      ji > 1 .AND. ji < jpi   &
252                      .AND. jj > 1 .AND. jj < jpj ) THEN
253                      zfwfad(jc) = zfwfad(jc) + emps_ad(ji,jj) /   &
254                          (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj))
255                      zfwfad(jc) = zfwfad(jc) + emp_ad(ji,jj)  /   &
256                          (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj))
257                  END IF
258                END DO
259            ELSE
260                DO jj = 1, jpj
261                   DO ji = 1, jpi
262                      zfwfad(jc) = zfwfad(jc) + emps_ad(ji,jj) / surf(jpncs+1)
263                      zfwfad(jc) = zfwfad(jc) + emp_ad(ji,jj)  / surf(jpncs+1)
264                   END DO
265                END DO
266            ENDIF
267         ELSEIF( ncstt(jc) == 2 ) THEN
268            ! Excess e-p+r (either sign) goes to open ocean, at outflow location
269            IF(      ji > 1 .AND. ji < jpi    &
270               .AND. jj > 1 .AND. jj < jpj ) THEN
271                DO jn = 1, ncsnr(jc)
272                  ji = mi0(ncsir(jc,jn))
273                  jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
274                  zfwfad(jc) = zfwfad(jc) + emps_ad(ji,jj)  &
275                      / (FLOAT(ncsnr(jc)) *  e1t(ji,jj) * e2t(ji,jj) )
276                  zfwfad(jc) = zfwfad(jc) + emp_ad(ji,jj)   &
277                      / (FLOAT(ncsnr(jc)) *  e1t(ji,jj) * e2t(ji,jj) )
278                END DO
279            ENDIF
280         ENDIF
281         !
282      END DO
283      !                                                   !--------------------!
284      !                                                   !  update emp, emps  !
285                                                          !--------------------!
286      IF( lk_mpp )   CALL mpp_sum ( zfwfad(:) , jpncs )       ! mpp: sum over all the global domain
287      DO jc = 1, jpncs
288         DO jj = ncsj1(jc), ncsj2(jc)
289            DO ji = ncsi1(jc), ncsi2(jc)
290               emp_ad(ji,jj) = emp_ad(ji,jj) + e1t(ji,jj) * e2t(ji,jj) * zfwfad(jc) * tmask_i(ji,jj)
291            END DO
292         END DO
293      END DO
294      zfwfad = 0.0_wp
295   END SUBROUTINE sbc_clo_adj
296   SUBROUTINE sbc_clo_adj_tst ( kumadt )
297      !!-----------------------------------------------------------------------
298      !!
299      !!                  ***  ROUTINE sbc_clo_adj_tst ***
300      !!
301      !! ** Purpose : Test the adjoint routine.
302      !!
303      !! ** Method  : Verify the scalar product
304      !!
305      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
306      !!
307      !!              where  L   = tangent routine
308      !!                     L^T = adjoint routine
309      !!                     W   = diagonal matrix of scale factors
310      !!                     dx  = input perturbation (random field)
311      !!                     dy  = L dx
312      !!
313      !! History :
314      !!        ! 09-08 (F. Vigilant)
315      !!-----------------------------------------------------------------------
316      !! * Modules used
317
318      !! * Arguments
319      INTEGER, INTENT(IN) :: &
320         & kumadt             ! Output unit
321
322      !! * Local declarations
323      INTEGER ::  &
324         & istp,  &
325         & jstp,  &
326         & ji,    &        ! dummy loop indices
327         & jj,    &
328         & jk
329      INTEGER, DIMENSION(jpi,jpj) :: &
330         & iseed_2d        ! 2D seed for the random number generator
331      REAL(KIND=wp) :: &
332         & zsp1,         & ! scalar product involving the tangent routine
333         & zsp2            ! scalar product involving the adjoint routine
334      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
335         & zemp_tlin ,     & ! Tangent input
336         & zemps_tlin ,    & ! Tangent input
337         & zemp_tlout,     & ! Tangent output
338         & zemps_tlout,    & ! Tangent output
339         & zemp_adin ,     & ! Adjoint input
340         & zemps_adin ,    & ! Adjoint input
341         & zemp_adout,     & ! Adjoint output
342         & zemps_adout,    & ! Adjoint output
343         & zr                ! 2D random field
344      CHARACTER(LEN=14) :: cl_name
345      ! Allocate memory
346
347      ALLOCATE( &
348         & zemp_tlin(  jpi,jpj),     &
349         & zemps_tlin( jpi,jpj),     &
350         & zemp_tlout( jpi,jpj),     &
351         & zemps_tlout(jpi,jpj),     &
352         & zemp_adin(  jpi,jpj),     &
353         & zemps_adin( jpi,jpj),     &
354         & zemp_adout( jpi,jpj),     &
355         & zemps_adout(jpi,jpj),     &
356         & zr(         jpi,jpj)      &
357         & )
358      !==================================================================
359      ! 1) dx = ( emp_tl, emps_tl, ssh_tl ) and
360      !    dy = ( emp_tl, emps_tl )
361      !==================================================================
362
363      !--------------------------------------------------------------------
364      ! Reset the tangent and adjoint variables
365      !--------------------------------------------------------------------
366          zemp_tlin  (:,:) = 0.0_wp
367          zemps_tlin (:,:) = 0.0_wp
368          zemp_tlout (:,:) = 0.0_wp
369          zemps_tlout(:,:) = 0.0_wp
370          zemp_adin  (:,:) = 0.0_wp
371          zemps_adin (:,:) = 0.0_wp
372          zemp_adout (:,:) = 0.0_wp
373          zemps_adout(:,:) = 0.0_wp
374          zr(:,:)          = 0.0_wp
375      !--------------------------------------------------------------------
376      ! Initialize the tangent input with random noise: dx
377      !--------------------------------------------------------------------
378
379      CALL grid_random( zr, 'T', 0.0_wp, stdemp )
380      DO jj = nldj, nlej
381         DO ji = nldi, nlei
382            zemp_tlin(ji,jj) = zr(ji,jj)
383         END DO
384      END DO
385      CALL grid_random(  zr, 'T', 0.0_wp, stdemp )
386      DO jj = nldj, nlej
387         DO ji = nldi, nlei
388            zemps_tlin(ji,jj) = zr(ji,jj)
389         END DO
390      END DO
391
392      emps_tl(:,:) = zemps_tlin(:,:)
393      emp_tl (:,:) = zemp_tlin (:,:)
394
395      CALL sbc_clo_tan( nit000 )
396
397      zemps_tlout(:,:) = emps_tl(:,:)
398      zemp_tlout (:,:) = emp_tl (:,:)
399      !-----------------------------------------------------------------
400      ! Initialize the adjoint variables: dy^* = W dy
401      !-----------------------------------------------------------------
402
403      DO jj = 1, jpi
404         DO ji = 1, jpj
405            zemp_adin( ji,jj) = zemp_tlout( ji,jj) &
406               &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) &
407               &               * tmask(ji,jj,1)
408            zemps_adin(ji,jj) = zemps_tlout(ji,jj) &
409               &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) &
410               &               * tmask(ji,jj,1)
411         END DO
412      END DO
413      !-----------------------------------------------------------------
414      ! Compute the scalar product: ( L dx )^T W dy
415      !-----------------------------------------------------------------
416      zsp1 = DOT_PRODUCT( zemp_tlout,  zemp_adin  ) &
417         & + DOT_PRODUCT( zemps_tlout, zemps_adin )
418
419      !-----------------------------------------------------------------
420      ! Call the adjoint routine: dx^* = L^T dy^*
421      !-----------------------------------------------------------------
422
423      emp_ad (:,:) = zemp_adin (:,:)
424      emps_ad(:,:) = zemps_adin(:,:)
425
426      CALL sbc_clo_adj ( nitend )
427
428      zemps_adout(:,:) = emps_ad(:,:)
429      zemp_adout (:,:) = emp_ad (:,:)
430
431      zsp2 = DOT_PRODUCT( zemp_tlin,  zemp_adout  ) &
432         & + DOT_PRODUCT( zemps_tlin, zemps_adout )
433
434      ! 14 char:'12345678901234'
435      cl_name = 'sbc_clo_adj   '
436      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
437
438      DEALLOCATE(       &
439         & zemp_tlin,   &
440         & zemps_tlin,  &
441         & zemp_tlout,  &
442         & zemps_tlout, &
443         & zemp_adin,   &
444         & zemps_adin,  &
445         & zemp_adout,  &
446         & zemps_adout, &
447         & zr           &
448         & )
449   END SUBROUTINE sbc_clo_adj_tst
450   !!======================================================================
451#endif
452END MODULE closea_tam
Note: See TracBrowser for help on using the repository browser.