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/TAM_V3_0/NEMOTAM/OPATAM_SRC/DOM – NEMO

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/DOM/closea_tam.F90 @ 2587

Last change on this file since 2587 was 1885, checked in by rblod, 14 years ago

add TAM sources

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