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

Last change on this file since 11931 was 11931, checked in by mathiot, 12 months ago

ENHANCE-02_ISF_nemo: add comments, improve memory usage of ln_isfcpl_cons option, fix issue in ISOMIP+ configuration

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