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

source: branches/TAM_V3_2_2/NEMOTAM/OPATAM_SRC/SBC/sbcssr_tam.F90 @ 3032

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

first import of NEMOTAM 3.2.2

File size: 16.4 KB
Line 
1MODULE sbcssr_tam
2#ifdef key_tam
3   !!======================================================================
4   !!                       ***  MODULE  sbcssr_tam  ***
5   !! Surface module :  add to heat and fresh water fluxes a restoring term
6   !!                   toward observed SST/SSS
7   !!                   Tangent and Adjoint Module
8   !!======================================================================
9   !! History of the direct routine:
10   !!            3.0  !  2006-06  (G. Madec)   Original code
11   !!            3.2  !  2009-04  (B. Lemaire) Introduce iom_put
12   !! History of the T&A routine:
13   !!            3.0  !  2008-11  (A. Vidard)  Original code (simplification: no linear salinity damping)
14   !!            3.2  !  2010-04  (A. Vidard)  Nemo3.2 update 
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   sbc_ssr        : add to sbc a restoring term toward SST/SSS climatology
19   !!----------------------------------------------------------------------
20   USE par_oce       , ONLY: & ! Ocean space and time domain variables
21      & jpi,                 &
22      & jpj,                 &
23      & jpiglo
24   USE par_kind      , ONLY: & ! Precision variables
25      & wp
26   USE dom_oce       , ONLY: & ! Ocean space and time domain
27      & e1t,                 &
28      & e2t,                 &
29# if defined key_vvl
30      & e3t_1,               &
31# else
32#  if defined key_zco
33      & e3t_0,               &
34#  else
35      & e3t,                 &
36#  endif
37# endif
38      & tmask,               &
39      & mig,                 &
40      & mjg,                 &
41      & nldi,                &
42      & nldj,                &
43      & nlei,                &
44      & nlej   
45   USE sbc_oce       , ONLY: & ! Surface boundary condition: ocean fields
46      & nn_fsbc
47   USE sbc_oce_tam   , ONLY: & ! Surface boundary condition: ocean fields
48      & sst_m_tl,            & ! mean (nn_fsbc time-step) surface sea temperature     [Celsius]
49      & sst_m_ad,            & ! mean (nn_fsbc time-step) surface sea temperature     [Celsius]
50      & qns_tl ,             & ! sea heat flux: non solar                 [W/m2]
51      & qns_ad             
52   USE in_out_manager, ONLY: & ! I/O manager
53      & lwp,                 &
54      & numout,              & 
55      & numnam,              &
56      & nit000,              & 
57      & nitend 
58   USE gridrandom    , ONLY: &     ! Random Gaussian noise on grids
59      & grid_random
60   USE dotprodfld,     ONLY: & ! Computes dot product for 3D and 2D fields
61      & dot_product
62   USE tstool_tam    , ONLY: &
63      & stdt,                &
64      & stdqns,              &
65      & prntst_adj
66
67   IMPLICIT NONE
68   PRIVATE
69
70   PUBLIC   sbc_ssr_tan    ! routine called in sbcmod_tam
71   PUBLIC   sbc_ssr_adj    ! routine called in sbcmod_tam
72   PUBLIC   sbc_ssr_adj_tst! routine called in tst
73
74   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  &
75     & erp_tl,    &    !: evaporation damping                          [kg/m2/s]
76     & qrp_tl,    &    !: heat flux damping                            [w/m2]
77     & erp_ad,    &    !: evaporation damping                          [kg/m2/s]
78     & qrp_ad          !: heat flux damping                            [w/m2]
79
80   !                                           !!* Namelist namsbc_ssr *
81   INTEGER, PUBLIC ::   nn_sstr     =   0       ! SST/SSS restoring indicator
82   INTEGER, PUBLIC ::   nn_sssr     =   0       ! SST/SSS restoring indicator
83   REAL(wp)        ::   rn_dqdt     = -40.e0    ! restoring factor on SST and SSS
84   REAL(wp)        ::   rn_deds     = -27.70    ! restoring factor on SST and SSS
85   LOGICAL         ::   ln_sssr_bnd = .false.   ! flag to bound erp term
86   REAL(wp)        ::   rn_sssr_bnd =   0.e0    ! ABS(Max./Min.) value of erp term [mm/day]
87
88   !! * Substitutions
89#  include "domzgr_substitute.h90"
90   !!----------------------------------------------------------------------
91   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
92   !! $Id$
93   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
94   !!----------------------------------------------------------------------
95
96CONTAINS
97
98   SUBROUTINE sbc_ssr_tan( kt )
99      !!---------------------------------------------------------------------
100      !!                     ***  ROUTINE sbc_ssr_tan  ***
101      !!
102      !! ** Purpose of the direct routine:
103      !!                Add to heat and/or freshwater fluxes a damping term
104      !!                toward observed SST and/or SSS.
105      !!
106      !! ** Method of the direct routine:
107      !!            : - Read namelist namsbc_ssr
108      !!              - Read observed SST and/or SSS
109      !!              - at each nscb time step
110      !!                   add a retroaction term on qns    (nn_sstr = 1)
111      !!                   add a damping term on emps       (nn_sssr = 1)
112      !!                   add a damping term on emp & emps (nn_sssr = 2)
113      !!---------------------------------------------------------------------
114      INTEGER, INTENT(in   ) ::   kt   ! ocean time step
115      !!
116      INTEGER  ::   ji, jj   ! dummy loop indices
117      REAL(wp) ::   zqrptl   ! local scalar for heat flux damping
118      !!
119      !!----------------------------------------------------------------------
120      !                                               ! -------------------- !
121      IF( kt == nit000 ) THEN                         ! First call kt=nit000 !
122         !                                            ! -------------------- !
123         CALL sbc_ssr_ini_tam ( 0 )
124      ENDIF
125
126      IF( nn_sstr + nn_sssr /= 0 ) THEN
127 
128         !                                         ! ========================= !
129         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN      !    Add restoring term     !
130            !                                      ! ========================= !
131            !
132            IF( nn_sstr == 1 ) THEN                   !* Temperature restoring term
133!CDIR COLLAPSE
134               DO jj = 1, jpj
135                  DO ji = 1, jpi
136                     zqrptl        = rn_dqdt * sst_m_tl(ji,jj)
137                     qns_tl(ji,jj) = qns_tl(ji,jj) + zqrptl
138                     qrp_tl(ji,jj) = zqrptl
139                  END DO
140               END DO
141            ENDIF
142            !
143            ! No linear Salinity damping term  (simplification)
144            !
145         ENDIF
146         !
147      ENDIF
148      !
149   END SUBROUTINE sbc_ssr_tan
150   SUBROUTINE sbc_ssr_adj( kt )
151      !!---------------------------------------------------------------------
152      !!                     ***  ROUTINE sbc_ssr_adj  ***
153      !!
154      !! ** Purpose of the direct routine:
155      !!                Add to heat and/or freshwater fluxes a damping term
156      !!                toward observed SST and/or SSS.
157      !!
158      !! ** Method of the direct routine:
159      !!            : - Read namelist namsbc_ssr
160      !!              - Read observed SST and/or SSS
161      !!              - at each nscb time step
162      !!                   add a retroaction term on qns    (nn_sstr = 1)
163      !!                   add a damping term on emps       (nn_sssr = 1)
164      !!                   add a damping term on emp & emps (nn_sssr = 2)
165      !!---------------------------------------------------------------------
166      INTEGER, INTENT(in   ) ::   kt   ! ocean time step
167      !!
168      INTEGER  ::   ji, jj   ! dummy loop indices
169      REAL(wp) ::   zqrpad   ! local scalar for heat flux damping
170      !!
171      !!----------------------------------------------------------------------
172      zqrpad = 0.0
173      !                                               ! -------------------- !
174      IF( kt == nitend ) THEN                         ! First call kt=nit000 !
175         !                                            ! -------------------- !
176         CALL sbc_ssr_ini_tam ( 1 )
177      ENDIF
178
179      IF( nn_sstr + nn_sssr /= 0 ) THEN
180 
181         !                                         ! ========================= !
182         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN      !    Add restoring term     !
183            !                                      ! ========================= !
184            !
185            IF( nn_sstr == 1 ) THEN                   ! Temperature restoring term
186!CDIR COLLAPSE
187               ! use zqrp scalar to optimize memory access (speedup the loop)
188               DO jj = 1, jpj
189                  DO ji = 1, jpi
190                     zqrpad = qrp_ad(ji,jj)
191                     qrp_ad(ji,jj) = 0.0_wp
192
193                     zqrpad = zqrpad + qns_ad(ji,jj)
194                     sst_m_ad(ji,jj) = sst_m_ad(ji,jj) + rn_dqdt * zqrpad
195                  END DO
196               END DO
197            ENDIF
198            !
199            ! No linear Salinity damping term  (simplification)
200            !
201         ENDIF
202         !
203      ENDIF
204      !
205   END SUBROUTINE sbc_ssr_adj
206   SUBROUTINE sbc_ssr_adj_tst( kumadt )
207      !!-----------------------------------------------------------------------
208      !!
209      !!                  ***  ROUTINE sbc_ssr_adj_tst ***
210      !!
211      !! ** Purpose : Test the adjoint routine.
212      !!
213      !! ** Method  : Verify the scalar product
214      !!           
215      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
216      !!
217      !!              where  L   = tangent routine
218      !!                     L^T = adjoint routine
219      !!                     W   = diagonal matrix of scale factors
220      !!                     dx  = input perturbation (random field)
221      !!                     dy  = L dx
222      !!
223      !!                   
224      !! History :
225      !!        ! 08-11 (A. Vidard)
226      !!-----------------------------------------------------------------------
227      !! * Modules used
228
229      !! * Arguments
230      INTEGER, INTENT(IN) :: &
231         & kumadt             ! Output unit
232 
233      !! * Local declarations
234      INTEGER ::  &
235         & ji,    &        ! dummy loop indices
236         & jj,    &       
237         & jk     
238      INTEGER, DIMENSION(jpi,jpj) :: &
239         & iseed_2d        ! 2D seed for the random number generator
240      REAL(KIND=wp) :: &
241         & zsp1,         & ! scalar product involving the tangent routine
242         & zsp2            ! scalar product involving the adjoint routine
243      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
244         & zsst_m_tlin ,   & ! Tangent input
245         & zqns_tlin ,     & ! Tangent input
246         & zqns_tlout,     & ! Tangent output
247         & zqrp_tlout,     & ! Tangent output
248         & zqns_adin ,     & ! Adjoint input
249         & zqrp_adin ,     & ! Adjoint input
250         & zsst_m_adout,   & ! Adjoint output
251         & zqns_adout,     & ! Adjoint output
252         & zr                ! 3D random field
253      CHARACTER(LEN=14) :: cl_name
254      ! Allocate memory
255
256      ALLOCATE( &
257         & zqns_tlin(   jpi,jpj),   &
258         & zsst_m_tlin( jpi,jpj),   &
259         & zqns_tlout(  jpi,jpj),   &
260         & zqrp_tlout(  jpi,jpj),   &
261         & zqns_adin(   jpi,jpj),   &
262         & zqrp_adin(   jpi,jpj),   &
263         & zqns_adout(  jpi,jpj),   &
264         & zsst_m_adout(jpi,jpj),   &
265         & zr(          jpi,jpj)    &
266         & )
267      !==================================================================
268      ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
269      !    dy = ( hdivb_tl, hdivn_tl )
270      !==================================================================
271
272      !--------------------------------------------------------------------
273      ! Reset the tangent and adjoint variables
274      !--------------------------------------------------------------------
275          zqns_tlin(   :,:) = 0.0_wp
276          zsst_m_tlin( :,:) = 0.0_wp
277          zqns_tlout(  :,:) = 0.0_wp
278          zqrp_tlout(  :,:) = 0.0_wp
279          zqns_adin(   :,:) = 0.0_wp
280          zqrp_adin(   :,:) = 0.0_wp
281          zqns_adout(  :,:) = 0.0_wp
282          zsst_m_adout(:,:) = 0.0_wp
283          zr(          :,:) = 0.0_wp
284      !--------------------------------------------------------------------
285      ! Initialize the tangent input with random noise: dx
286      !--------------------------------------------------------------------
287
288      DO jj = 1, jpj
289         DO ji = 1, jpi
290            iseed_2d(ji,jj) = - ( 596035 + &
291               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
292         END DO
293      END DO
294      CALL grid_random( iseed_2d, zr, 'T', 0.0_wp, stdqns )
295      DO jj = nldj, nlej
296         DO ji = nldi, nlei
297            zqns_tlin(ji,jj) = zr(ji,jj) 
298         END DO
299      END DO
300
301      DO jj = 1, jpj
302         DO ji = 1, jpi
303            iseed_2d(ji,jj) = - ( 371498 + &
304               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
305         END DO
306      END DO
307      CALL grid_random( iseed_2d, zr, 'T', 0.0_wp, stdt )
308      DO jj = nldj, nlej
309         DO ji = nldi, nlei
310            zsst_m_tlin(ji,jj) = zr(ji,jj) 
311         END DO
312      END DO
313
314      sst_m_tl(:,:) = zsst_m_tlin(:,:)
315      qns_tl(  :,:) = zqns_tlin(  :,:)
316
317      CALL sbc_ssr_tan (nit000)
318      zqns_tlout(:,:) = qns_tl(:,:)
319      zqrp_tlout(:,:) = qrp_tl(:,:)
320      !--------------------------------------------------------------------
321      ! Initialize the adjoint variables: dy^* = W dy
322      !--------------------------------------------------------------------
323
324        DO jj = nldj, nlej
325           DO ji = nldi, nlei
326              zqns_adin(ji,jj)   = zqns_tlout(ji,jj) &
327                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) &
328                 &               * tmask(ji,jj,1)
329            END DO
330         END DO
331        DO jj = nldj, nlej
332           DO ji = nldi, nlei
333              zqrp_adin(ji,jj)   = zqrp_tlout(ji,jj) &
334                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) &
335                 &               * tmask(ji,jj,1)
336            END DO
337         END DO
338      !--------------------------------------------------------------------
339      ! Compute the scalar product: ( L dx )^T W dy
340      !--------------------------------------------------------------------
341
342      zsp1 = DOT_PRODUCT( zqns_tlout, zqns_adin ) &
343         & + DOT_PRODUCT( zqrp_tlout, zqrp_adin ) 
344
345      !--------------------------------------------------------------------
346      ! Call the adjoint routine: dx^* = L^T dy^*
347      !--------------------------------------------------------------------
348      qns_ad(:,:) = zqns_adin(:,:)
349      qrp_ad(:,:) = zqrp_adin(:,:)
350      CALL sbc_ssr_adj (nit000)
351      zqns_adout(  :,:) = qns_ad(  :,:)
352      zsst_m_adout(:,:) = sst_m_ad(:,:)
353
354      zsp2 = DOT_PRODUCT( zqns_tlin,   zqns_adout   ) &
355         & + DOT_PRODUCT( zsst_m_tlin, zsst_m_adout )
356
357      ! 14 char:'12345678901234'
358      cl_name = 'sbc_ssr_adj   '
359      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
360
361      DEALLOCATE(        &
362         & zqns_tlin,    &
363         & zsst_m_tlin,  &
364         & zqns_tlout,   &
365         & zqrp_tlout,   &
366         & zqns_adin,    &
367         & zqrp_adin,    &
368         & zqns_adout,   &
369         & zsst_m_adout, &
370         & zr            &
371         & )
372
373
374
375   END SUBROUTINE sbc_ssr_adj_tst
376
377
378   SUBROUTINE sbc_ssr_ini_tam( kindic )
379      USE fldread
380      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
381      TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read
382      !!----------------------------------------------------------------------
383      INTEGER, INTENT(IN) :: kindic
384      NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, &
385         &                 sn_sss, ln_sssr_bnd, rn_sssr_bnd
386
387      REWIND ( numnam )         ! ... read in namlist namflx
388      READ( numnam, namsbc_ssr ) 
389
390      IF(lwp) THEN              ! control print
391         WRITE(numout,*)
392         WRITE(numout,*) 'sbc_ssr_tam : SST and/or SSS damping term '
393         WRITE(numout,*) '~~~~~~~~~~~ '
394         WRITE(numout,*) '          SST restoring term (Yes=1)             nn_sstr = ', nn_sstr
395         WRITE(numout,*) '          SSS damping term (Yes=1, salt flux)    nn_sssr = ', nn_sssr
396         WRITE(numout,*) '                           (Yes=2, volume flux) '
397         WRITE(numout,*) '          dQ/dT (restoring magnitude on SST)     dqdt    = ', rn_dqdt, ' W/m2/K'
398         WRITE(numout,*) '          dE/dS (restoring magnitude on SST)     deds    = ', rn_deds, ' mm/day'
399      ENDIF
400     
401      !
402      ! Initialize qrp and erp if no restoring
403      IF ( kindic == 0 ) THEN
404         qrp_tl(:,:) = 0.e0 
405         erp_tl(:,:) = 0.e0
406      ELSEIF ( kindic == 1 ) THEN
407         qrp_ad(:,:) = 0.e0 
408         erp_ad(:,:) = 0.e0 
409      END IF
410   END SUBROUTINE sbc_ssr_ini_tam
411
412   !!======================================================================
413#endif
414END MODULE sbcssr_tam
Note: See TracBrowser for help on using the repository browser.