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/dev_r11943_MERGE_2019/src/OCE/ISF – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isftbl.F90 @ 12242

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

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