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.
isftbl.F90 in NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF – NEMO

source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isftbl.F90 @ 11844

Last change on this file since 11844 was 11844, checked in by mathiot, 4 years ago

rm useless lbclnk + fix reproducibility WED025 + add isf debug option
File size: 12.9 KB
Line 
1MODULE isftbl
2   !!======================================================================
3   !!                       ***  MODULE  isftbl  ***
4   !! isftbl module :  compute properties of top boundary layer
5   !!======================================================================
6   !! History :  4.1  !  2019-09  (P. Mathiot) original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   isftbl       : routine to compute :
11   !!                  - geometry of the ice shelf tbl (isf_tbl_lvl, isftbl_ktop, isftbl_kbot)
12   !!                    (top and bottom level, thickness and fraction of deepest level affected)
13   !!                  - tbl averaged properties (isf_tbl, isf_tbl_avg)
14   !!----------------------------------------------------------------------
15
16   USE dom_oce ! vertical scale factor
17   USE lbclnk  ! lbc_lnk subroutine
18   USE isfutils
19   USE isf
20
21   IMPLICIT NONE
22
23   PRIVATE
24
25   PUBLIC isf_tbl, isf_tbl_avg, isf_tbl_lvl, isf_tbl_ktop, isf_tbl_kbot
26
27CONTAINS
28
29   SUBROUTINE isf_tbl( pvarin, pvarout, cd_ptin, ktop, phtbl, kbot, pfrac )
30      !!--------------------------------------------------------------------
31      !!                  ***  SUBROUTINE isf_tbl  ***
32      !!
33      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
34      !!
35      !! ** Method : Average properties over a specific thickness
36      !!
37      !! ** Reference : inspired from : Losch, Modeling ice shelf cavities in a z coordinate ocean general circulation model
38      !!                https://doi.org/10.1029/2007JC004368 , 2008
39      !!
40      !!--------------------------------------------------------------------
41      !!-------------------------- OUT -------------------------------------
42      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(  out) :: pvarout ! 2d average of pvarin
43      !!-------------------------- IN  -------------------------------------
44      CHARACTER(len=1)                      , INTENT(in   ) :: cd_ptin       ! point of variable in/out
45      REAL(wp), DIMENSION(jpi,jpj,jpk)      , INTENT(in   ) :: pvarin        ! 3d variable to average over the tbl
46      INTEGER,  DIMENSION(jpi,jpj)          , INTENT(in   ) :: ktop          ! top level
47      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in   ) :: phtbl         ! tbl thickness
48      !!-------------------------- IN OPTIONAL -----------------------------
49      INTEGER,  DIMENSION(jpi,jpj), OPTIONAL, INTENT(in   ) :: kbot          ! bottom level
50      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in   ) :: pfrac         ! fraction of bottom cell affected by tbl
51      !!--------------------------------------------------------------------
52      INTEGER ::   ji, jj                     ! loop index
53      INTEGER , DIMENSION(jpi,jpj) :: ikbot   ! bottom level of the tbl
54      REAL(wp), DIMENSION(jpi,jpj) :: zvarout ! 2d average of pvarin
55      REAL(wp), DIMENSION(jpi,jpj) :: zhtbl   ! thickness of the tbl
56      REAL(wp), DIMENSION(jpi,jpj) :: zfrac   ! thickness of the tbl
57      !!--------------------------------------------------------------------
58      !
59      SELECT CASE ( cd_ptin )
60      CASE ( 'U' )
61         !
62         ! copy phtbl (phtbl is INTENT in as we don't want to change it)
63         zhtbl = phtbl
64         !
65         ! compute tbl lvl and thickness
66         CALL isf_tbl_lvl( hu_n, e3u_n, ktop, ikbot, zhtbl, zfrac )
67         !
68         ! compute tbl property at U point
69         CALL isf_tbl_avg( miku, ikbot, zhtbl, zfrac, e3u_n, pvarin, zvarout )
70         !
71         ! check if needed (probably yes)
72         CALL lbc_lnk('sbcisf', pvarout,'U',-1.)
73         !
74         ! compute tbl property at T point
75         pvarout(1,:) = 0._wp
76         DO jj = 1, jpj
77            DO ji = 2, jpi
78               pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji-1,jj))
79            END DO
80         END DO
81         ! lbclnk not needed as a final communication is done after the computation of fwf
82         !
83      CASE ( 'V' )
84         !
85         ! copy phtbl (phtbl is INTENT in as we don't want to change it)
86         zhtbl = phtbl
87         !
88         ! compute tbl lvl and thickness
89         CALL isf_tbl_lvl( hv_n, e3v_n, ktop, ikbot, zhtbl, zfrac )
90         !
91         ! compute tbl property at V point
92         CALL isf_tbl_avg( mikv, ikbot, zhtbl, zfrac, e3v_n, pvarin, zvarout )
93         !
94         ! check if needed (probably yes)
95         CALL lbc_lnk('sbcisf', pvarout,'V',-1.)
96         !
97         ! pvarout is an averaging of wet point
98         pvarout(:,1) = 0._wp
99         DO jj = 2, jpj
100            DO ji = 1, jpi
101               pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji,jj-1))
102            END DO
103         END DO
104         ! lbclnk not needed as a final communication is done after the computation of fwf
105         !
106      CASE ( 'T' )
107         !
108         ! compute tbl property at T point
109         CALL isf_tbl_avg( ktop, kbot, phtbl, pfrac, e3t_n, pvarin, pvarout )
110         !
111      END SELECT
112      !
113      IF (ln_isfdebug) CALL debug('isf_tbl pvarout:',pvarout)
114      !
115   END SUBROUTINE isf_tbl
116
117   SUBROUTINE isf_tbl_avg( ktop, kbot, phtbl, pfrac, pe3, pvarin, pvarout )
118      !!--------------------------------------------------------------------
119      !!                  ***  ROUTINE isf_tbl_avg  ***
120      !!
121      !! ** Purpose : compute mean property in the boundary layer
122      !!
123      !! ** Method  : Depth average is made between the top level ktop and the bottom level kbot
124      !!              over a thickness phtbl. The bottom level is partially counted (pfrac).
125      !!
126      !!--------------------------------------------------------------------
127      !!-------------------------- OUT -------------------------------------
128      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out) :: pvarout      ! tbl property averaged over phtbl between level ktop and kbot
129      !!-------------------------- IN  -------------------------------------
130      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(in   ) :: ktop, kbot   ! top and bottom level of the top boundary layer
131      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: phtbl, pfrac ! fraction of bottom level to be affected by the tbl
132      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pe3          ! vertical scale factor
133      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pvarin       ! tbl property to average between ktop, kbot over phtbl
134      !!--------------------------------------------------------------------
135      INTEGER  :: ji,jj,jk                    ! loop indices
136      INTEGER  :: ikt, ikb                    ! top and bottom levels
137      !!--------------------------------------------------------------------
138      !
139      ! compute tbl top.bottom level and thickness
140      DO jj = 1,jpj
141         DO ji = 1,jpi
142            !
143            ! tbl top/bottom indices initialisation
144            ikt = ktop(ji,jj) ; ikb = kbot(ji,jj)
145            !
146            ! level fully include in the ice shelf boundary layer
147            pvarout(ji,jj) = SUM( pvarin(ji,jj,ikt:ikb-1) * pe3(ji,jj,ikt:ikb-1) ) / phtbl(ji,jj)
148            !
149            ! level partially include in ice shelf boundary layer
150            pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * pe3(ji,jj,ikb) / phtbl(ji,jj) * pfrac(ji,jj)
151            !
152         END DO
153      END DO
154
155   END SUBROUTINE isf_tbl_avg
156
157   SUBROUTINE isf_tbl_lvl( phw, pe3, ktop, kbot, phtbl, pfrac )
158      !!--------------------------------------------------------------------
159      !!                  ***  ROUTINE isf_tbl_lvl  ***
160      !!
161      !! ** Purpose : - compute bottom level off the top boundary layer
162      !!              - thickness of the top boundary layer
163      !!              - fraction of the bottom level affected by the tbl
164      !!
165      !!---------------------------------------------------------------------
166      !!-------------------------- OUT --------------------------------------
167      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(  out) :: kbot   ! bottom level of the top boundary layer
168      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out) :: pfrac  ! fraction of bottom level in the tbl
169      !!-------------------------- IN  --------------------------------------
170      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(in   ) :: ktop   ! top level of the top boundary layer
171      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: phw    ! water column thickness
172      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pe3    ! vertical scale factor
173      !!-------------------------- INOUT ------------------------------------
174      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(inout) :: phtbl  ! top boundary layer thickness
175      !!---------------------------------------------------------------------
176      INTEGER :: ji,jj,jk
177      INTEGER :: ikt, ikb
178      !!---------------------------------------------------------------------
179      !
180      ! get htbl
181      DO jj = 1,jpj
182         DO ji = 1,jpi
183            !
184            ! tbl top/bottom indices initialisation
185            ikt = ktop(ji,jj)
186            !
187            ! limit the tbl to water thickness.
188            phtbl(ji,jj) = MIN( phtbl(ji,jj), phw(ji,jj) )
189            !
190            ! thickness of boundary layer must be at least the top level thickness
191            phtbl(ji,jj) = MAX( phtbl(ji,jj), pe3(ji,jj,ikt) )
192            !
193         END DO
194      END DO
195      !
196      ! get ktbl
197      CALL isf_tbl_kbot(ktop, phtbl, pe3, kbot)
198      !
199      ! get pfrac
200      DO jj = 1,jpj
201         DO ji = 1,jpi
202            !
203            ! tbl top/bottom indices initialisation
204            ikt = ktop(ji,jj) ; ikb = kbot(ji,jj)
205            !
206            ! proportion of the bottom cell included in ice shelf boundary layer
207            pfrac(ji,jj) = ( phtbl(ji,jj) - SUM( pe3(ji,jj,ikt:ikb-1) ) ) / pe3(ji,jj,ikb)
208            !
209         END DO
210      END DO
211      !
212   END SUBROUTINE isf_tbl_lvl
213   !
214   SUBROUTINE isf_tbl_kbot(ktop, phtbl, pe3, kbot)
215      !!--------------------------------------------------------------------
216      !!                  ***  ROUTINE isf_tbl_bot  ***
217      !!
218      !! ** Purpose : compute bottom level of the isf top boundary layer
219      !!
220      !!--------------------------------------------------------------------
221      !!-------------------------- OUT -------------------------------------
222      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(  out) :: kbot   ! bottom level of the top boundary layer
223      !!-------------------------- IN  -------------------------------------
224      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: phtbl  ! top boundary layer thickness
225      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(in   ) :: ktop   ! top level of the top boundary layer
226      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pe3    ! vertical scale factor
227      !!--------------------------------------------------------------------
228      INTEGER :: ji, jj
229      INTEGER :: ikt, ikb
230      !!--------------------------------------------------------------------
231      !
232      ! phtbl need to be bounded by water column thickness before
233      ! test: if htbl = water column thickness, should return mbathy
234      ! test: if htbl = 0 should return ktop (phtbl cap to e3t(ji,jj,1))
235      !
236      ! get ktbl
237      DO jj = 1,jpj
238         DO ji = 1,jpi
239            !
240            ! determine the deepest level influenced by the boundary layer
241            ikt = ktop(ji,jj)
242            ikb = ikt
243            DO WHILE ( SUM(pe3(ji,jj,ikt:ikb-1)) < phtbl(ji,jj ) ) ;  ikb = ikb + 1 ;  END DO
244            kbot(ji,jj) = ikb - 1
245            !
246         END DO
247      END DO
248      !
249   END SUBROUTINE isf_tbl_kbot
250      !
251   SUBROUTINE isf_tbl_ktop(pdep, ktop)
252      !!--------------------------------------------------------------------
253      !!                  ***  ROUTINE isf_tbl_top  ***
254      !!
255      !! ** Purpose : compute top level of the isf top boundary layer in case of an ice shelf parametrisation
256      !!
257      !!--------------------------------------------------------------------
258      !!-------------------------- OUT -------------------------------------
259      INTEGER,  DIMENSION(jpi,jpj), INTENT(  out) :: ktop        ! top level affected by the ice shelf parametrisation
260      !!-------------------------- IN  -------------------------------------
261      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pdep        ! top depth of the parametrisation influence
262      !!--------------------------------------------------------------------
263      INTEGER :: ji,jj
264      INTEGER :: ikt
265      !!--------------------------------------------------------------------
266      !
267      ! compute top level (need to be recomputed each time (z*, z~)
268      ! be sure pdep is already correctly bounded
269      ! test: this routine run on isfdraft should return mikt
270      ! test: this routine run with pdep = 0 should return 1
271      !
272      DO ji = 1, jpi
273         DO jj = 1, jpj
274            ikt = 2
275            DO WHILE ( gdepw_n(ji,jj,ikt) <= pdep(ji,jj ) ) ;  ikt = ikt + 1 ;  END DO
276            ktop(ji,jj) = ikt - 1
277         END DO
278      END DO
279      !
280   END SUBROUTINE isf_tbl_ktop
281
282END MODULE isftbl
Note: See TracBrowser for help on using the repository browser.