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_0/NEMOTAM/OPATAM_SRC/SBC – NEMO

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/SBC/sbcssr_tam.F90 @ 1885

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

add TAM sources

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