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_YRR_v2.F90 on Ticket #2071 – Attachment – NEMO

Ticket #2071: sbcssr_YRR_v2.F90

File sbcssr_YRR_v2.F90, 15.2 KB (added by yruprich, 6 years ago)
Line 
1MODULE sbcssr
2   !!======================================================================
3   !!                       ***  MODULE  sbcssr  ***
4   !! Surface module :  heat and fresh water flux restoring term toward targeted SST/SSS
5   !!======================================================================
6   !! History :  3.0  !  2006-06  (G. Madec)  Original code
7   !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put
8   !!            3.6  !  2018-02  (Y. Ruprich-Robert) restoring ponderated by sea-ice fraction
9   !!            3.6  !  2018-02  (Y. Ruprich-Robert) Subregion restoring mask
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   sbc_ssr       : add to sbc a restoring term toward SST/SSS target
14   !!   sbc_ssr_init  : initialisation of surface restoring
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and tracers
17   USE dom_oce        ! ocean space and time domain
18   USE sbc_oce        ! surface boundary condition
19   USE phycst         ! physical constants
20   USE sbcrnf         ! surface boundary condition : runoffs
21   !
22   USE fldread        ! read input fields
23   USE iom            ! I/O manager
24   USE in_out_manager ! I/O manager
25   USE lib_mpp        ! distribued memory computing library
26   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
27   USE timing         ! Timing
28   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   sbc_ssr        ! routine called in sbcmod
34   PUBLIC   sbc_ssr_init   ! routine called in sbcmod
35
36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   erp   !: evaporation damping   [kg/m2/s]
37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qrp   !: heat flux damping        [w/m2]
38! 02/2018 - Yohan Ruprich-Robert change: subregion restoring mask
39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sf_msk_f        !: restoring mask
40   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sf_msk_f_init   !: initial restoring mask
41
42
43   !                                   !!* Namelist namsbc_ssr *
44   INTEGER, PUBLIC ::   nn_sstr         ! SST restoring indicator
45   INTEGER, PUBLIC ::   nn_sssr         ! SSS restoring indicator
46! 02/2018 - Yohan Ruprich-Robert change: no restoring where sea-ice
47   INTEGER, PUBLIC ::   nn_icer         ! SST/SSS restoring indicator where sea-ice
48! 02/2018 - Yohan Ruprich-Robert change: subregion restoring mask
49   INTEGER, PUBLIC ::   nn_msk          ! SST/SSS restoring mask indicator
50   REAL(wp)        ::   rn_dqdt         ! restoring factor on SST and SSS
51   REAL(wp)        ::   rn_deds         ! restoring factor on SST and SSS
52   LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term
53   REAL(wp)        ::   rn_sssr_bnd     ! ABS(Max./Min.) value of erp term [mm/day]
54
55   REAL(wp) , ALLOCATABLE, DIMENSION(:) ::   buffer   ! Temporary buffer for exchange
56   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sst   ! structure of input SST (file informations, fields read)
57   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sss   ! structure of input SSS (file informations, fields read)
58   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_msk   ! structure of input MASK (file informations, fields read)
59
60   !! * Substitutions
61#  include "domzgr_substitute.h90"
62   !!----------------------------------------------------------------------
63   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
64   !! $Id: sbcssr.F90 4990 2014-12-15 16:42:49Z timgraham $
65   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
66   !!----------------------------------------------------------------------
67CONTAINS
68
69   SUBROUTINE sbc_ssr( kt )
70      !!---------------------------------------------------------------------
71      !!                     ***  ROUTINE sbc_ssr  ***
72      !!
73      !! ** Purpose :   Add to heat and/or freshwater fluxes a damping term
74      !!                toward targeted SST and/or SSS.
75      !!
76      !! ** Method  : - Read namelist namsbc_ssr
77      !!              - Read targeted SST and/or SSS
78      !!              - at each nscb time step
79      !!                   add a retroaction term on qns    (nn_sstr = 1)
80      !!                   add a damping term on sfx        (nn_sssr = 1)
81      !!                   add a damping term on emp        (nn_sssr = 2)
82      !!---------------------------------------------------------------------
83      INTEGER, INTENT(in   ) ::   kt   ! ocean time step
84      !!
85      INTEGER  ::   ji, jj   ! dummy loop indices
86      REAL(wp) ::   zerp     ! local scalar for evaporation damping
87      REAL(wp) ::   zqrp     ! local scalar for heat flux damping
88      REAL(wp) ::   zsrp     ! local scalar for unit conversion of rn_deds factor
89      REAL(wp) ::   zerp_bnd ! local scalar for unit conversion of rn_epr_max factor
90      INTEGER  ::   ierror   ! return error code
91      !!
92      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
93      TYPE(FLD_N) ::   sn_sst, sn_sss, sn_msk    ! informations about the fields to be read
94      !!----------------------------------------------------------------------
95      !     
96      IF( nn_timing == 1 )  CALL timing_start('sbc_ssr')
97      !
98      IF( nn_sstr + nn_sssr /= 0 ) THEN
99         !
100! 02/2018 - Yohan Ruprich-Robert change: read mask and provides it at kt
101! To be activated only if the mask changes along the simulation (e.g., seasonal forcing)
102!         CALL fld_read( kt, nn_fsbc, sf_msk )   ! Read mask data and provides it at kt
103!         sf_msk_f_init = sf_msk(1)%fnow(:,:,1)
104
105! 02/2018 - Yohan Ruprich-Robert change: modify mask based on sea-ice restoring option
106         sf_msk_f = sf_msk_f_init      ! initialize mask at time step kt
107         IF( nn_icer == 0 ) THEN       ! case surface restoring ponderate by sea-ice fraction
108            sf_msk_f = sf_msk_f(:,:) * ( 1.e0 - fr_i(:,:) )
109         ENDIF
110         !
111         IF( nn_sstr == 1)   CALL fld_read( kt, nn_fsbc, sf_sst )   ! Read SST data and provides it at kt
112         IF( nn_sssr >= 1)   CALL fld_read( kt, nn_fsbc, sf_sss )   ! Read SSS data and provides it at kt
113         !
114         !                                         ! ========================= !
115         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN      !    Add restoring term     !
116            !                                      ! ========================= !
117           IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term
118             DO jj = 1, jpj
119                  DO ji = 1, jpi
120! 02/2018 - Yohan Ruprich-Robert change: subregion restoring mask
121                     zqrp = sf_msk_f(ji,jj) * rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) )
122                     qns(ji,jj) = qns(ji,jj) + zqrp
123                     qrp(ji,jj) = zqrp
124                  END DO
125               END DO
126            ENDIF
127            CALL iom_put( "qrp", qrp )
128
129            !
130            IF( nn_sssr == 1 ) THEN                                   !* Salinity damping term (salt flux only (sfx))
131               zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
132!CDIR COLLAPSE
133               DO jj = 1, jpj
134                  DO ji = 1, jpi
135! 02/2018 - Yohan Ruprich-Robert change: subregion restoring mask
136                     zerp = sf_msk_f(ji,jj) * zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
137                          &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) 
138                     sfx(ji,jj) = sfx(ji,jj) + zerp                  ! salt flux
139                     erp(ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic only)
140                  END DO
141               END DO
142               !
143            ELSEIF( nn_sssr == 2 ) THEN                               !* Salinity damping term (volume flux (emp) and associated heat flux (qns)
144               zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s]
145               zerp_bnd = rn_sssr_bnd / rday                          !       -              -   
146!CDIR COLLAPSE
147               DO jj = 1, jpj
148                  DO ji = 1, jpi                           
149! 02/2018 - Yohan Ruprich-Robert change: subregion restoring mask
150                     zerp = sf_msk_f(ji,jj) * zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
151                          &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   &
152                          &        / MAX(  sss_m(ji,jj), 1.e-20   )
153                     IF( ln_sssr_bnd )   zerp = SIGN( 1., zerp ) * MIN( zerp_bnd, ABS(zerp) )
154                     emp(ji,jj) = emp (ji,jj) + zerp
155                     qns(ji,jj) = qns(ji,jj) - zerp * rcp * sst_m(ji,jj)
156                     erp(ji,jj) = zerp
157                  END DO
158               END DO
159            ENDIF
160            CALL iom_put( "erp", erp )                             ! freshwater flux damping
161
162         ENDIF
163            !
164      ENDIF
165         !
166      IF( nn_timing == 1 )  CALL timing_stop('sbc_ssr')
167      !
168   END SUBROUTINE sbc_ssr
169
170 
171   SUBROUTINE sbc_ssr_init
172      !!---------------------------------------------------------------------
173      !!                  ***  ROUTINE sbc_ssr_init  ***
174      !!
175      !! ** Purpose :   initialisation of surface damping term
176      !!
177      !! ** Method  : - Read namelist namsbc_ssr
178      !!              - Read observed SST and/or SSS if required
179      !!---------------------------------------------------------------------
180      INTEGER  ::   ji, jj   ! dummy loop indices
181      REAL(wp) ::   zerp     ! local scalar for evaporation damping
182      REAL(wp) ::   zqrp     ! local scalar for heat flux damping
183      REAL(wp) ::   zsrp     ! local scalar for unit conversion of rn_deds factor
184      REAL(wp) ::   zerp_bnd ! local scalar for unit conversion of rn_epr_max factor
185      INTEGER  ::   ierror   ! return error code
186      LOGICAL  ::   llok
187      !!
188      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
189      TYPE(FLD_N) ::   sn_sst, sn_sss, sn_msk   ! informations about the fields to be read
190      NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, nn_icer, nn_msk, rn_dqdt, rn_deds, sn_sst, sn_sss, sn_msk, ln_sssr_bnd, rn_sssr_bnd
191      INTEGER     ::  ios
192      !!----------------------------------------------------------------------
193      !
194 
195      REWIND( numnam_ref )              ! Namelist namsbc_ssr in reference namelist :
196      READ  ( numnam_ref, namsbc_ssr, IOSTAT = ios, ERR = 901)
197901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_ssr in reference namelist', lwp )
198
199      REWIND( numnam_cfg )              ! Namelist namsbc_ssr in configuration namelist :
200      READ  ( numnam_cfg, namsbc_ssr, IOSTAT = ios, ERR = 902 )
201902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_ssr in configuration namelist', lwp )
202      IF(lwm) WRITE ( numond, namsbc_ssr )
203
204      IF(lwp) THEN                 !* control print
205         WRITE(numout,*)
206         WRITE(numout,*) 'sbc_ssr : SST and/or SSS damping term '
207         WRITE(numout,*) '~~~~~~~ '
208         WRITE(numout,*) '   Namelist namsbc_ssr :'
209         WRITE(numout,*) '      SST restoring term (Yes=1)             nn_sstr     = ', nn_sstr
210         WRITE(numout,*) '      SSS damping term (Yes=1, salt flux)    nn_sssr     = ', nn_sssr
211         WRITE(numout,*) '                       (Yes=2, volume flux) '
212         WRITE(numout,*) '      restoring where sea-ice (Yes=1)        nn_icer     = ', nn_icer
213         WRITE(numout,*) '      subregion restoring mask (Yes=1)       nn_msk      = ', nn_msk
214         WRITE(numout,*) '      dQ/dT (restoring magnitude on SST)     rn_dqdt     = ', rn_dqdt, ' W/m2/K'
215         WRITE(numout,*) '      dE/dS (restoring magnitude on SST)     rn_deds     = ', rn_deds, ' mm/day'
216         WRITE(numout,*) '      flag to bound erp term                 ln_sssr_bnd = ', ln_sssr_bnd
217         WRITE(numout,*) '      ABS(Max./Min.) erp threshold           rn_sssr_bnd = ', rn_sssr_bnd, ' mm/day'
218      ENDIF
219      !
220      !                            !* Allocate erp and qrp array
221      ALLOCATE( qrp(jpi,jpj), erp(jpi,jpj), sf_msk_f_init(jpi,jpj), sf_msk_f(jpi,jpj), STAT=ierror )
222      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate erp and qrp array' )
223      !
224! 02/2018 - Yohan Ruprich-Robert change: subregion restoring mask
225      IF( nn_sstr + nn_sssr /= 0 ) THEN
226         !
227         sf_msk_f_init = tmask(:,:,1)
228         IF (nn_msk == 1) THEN             !* set sf_msk structure & allocate arrays
229
230            ALLOCATE( sf_msk(1), STAT=ierror )
231            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_msk structure' )
232            ALLOCATE( sf_msk(1)%fnow(jpi,jpj,1), STAT=ierror )
233            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_msk now array' )
234            !
235            ! fill sf_msk with sn_msk and control print
236            CALL fld_fill( sf_msk, (/ sn_msk /), cn_dir, 'sbc_ssr', 'mask for sea surface restoring', 'namsbc_ssr' )
237            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_msk data array' )
238            !
239! 02/2018 - Yohan Ruprich-Robert change: read mask
240! Comment the following 2 lines if the mask changes along the simulation
241            CALL fld_read( 0, nn_fsbc, sf_msk )   ! Read mask
242            sf_msk_f_init = sf_msk(1)%fnow(:,:,1)
243         ENDIF
244      ENDIF
245
246      IF( nn_sstr == 1 ) THEN      !* set sf_sst structure & allocate arrays
247         !
248         ALLOCATE( sf_sst(1), STAT=ierror )
249         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst structure' )
250         ALLOCATE( sf_sst(1)%fnow(jpi,jpj,1), STAT=ierror )
251         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst now array' )
252         !
253         ! fill sf_sst with sn_sst and control print
254         CALL fld_fill( sf_sst, (/ sn_sst /), cn_dir, 'sbc_ssr', 'SST restoring term toward SST data', 'namsbc_ssr' )
255         IF( sf_sst(1)%ln_tint )   ALLOCATE( sf_sst(1)%fdta(jpi,jpj,1,2), STAT=ierror )
256         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst data array' )
257         !
258      ENDIF
259      !
260      IF( nn_sssr >= 1 ) THEN      !* set sf_sss structure & allocate arrays
261         !
262         ALLOCATE( sf_sss(1), STAT=ierror )
263         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss structure' )
264         ALLOCATE( sf_sss(1)%fnow(jpi,jpj,1), STAT=ierror )
265         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss now array' )
266         !
267         ! fill sf_sss with sn_sss and control print
268         CALL fld_fill( sf_sss, (/ sn_sss /), cn_dir, 'sbc_ssr', 'SSS restoring term toward SSS data', 'namsbc_ssr' )
269         IF( sf_sss(1)%ln_tint )   ALLOCATE( sf_sss(1)%fdta(jpi,jpj,1,2), STAT=ierror )
270         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss data array' )
271         !
272      ENDIF
273      !
274      !                            !* Initialize qrp and erp if no restoring
275      IF( nn_sstr /= 1                   )   qrp(:,:) = 0._wp
276      IF( nn_sssr /= 1 .OR. nn_sssr /= 2 )   erp(:,:) = 0._wp
277      !
278   END SUBROUTINE sbc_ssr_init
279     
280   !!======================================================================
281END MODULE sbcssr