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 @ 11541

Last change on this file since 11541 was 11541, checked in by mathiot, 5 years ago

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