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

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/ISF/isftbl.F90 @ 13151

Last change on this file since 13151 was 13151, checked in by gm, 4 years ago

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