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

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

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