source: NEMO/trunk/src/OCE/ISF/isfhdiv.F90 @ 12489

Last change on this file since 12489 was 12489, checked in by davestorkey, 7 months ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

File size: 6.4 KB
Line 
1MODULE isfhdiv
2   !!======================================================================
3   !!                       ***  MODULE  isfhdiv  ***
4   !! ice shelf horizontal divergence module :  update the horizontal divergence
5   !!                   with the ice shelf melt and coupling correction
6   !!======================================================================
7   !! History :  4.0  !  2019-09  (P. Mathiot) Original code
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   isf_hdiv    : update the horizontal divergence with the ice shelf
12   !!                 melt and coupling correction
13   !!----------------------------------------------------------------------
14
15   USE isf_oce                ! ice shelf
16
17   USE dom_oce                ! time and space domain
18   USE phycst , ONLY: r1_rho0 ! physical constant
19   USE in_out_manager         !
20
21   IMPLICIT NONE
22
23   PRIVATE
24
25   PUBLIC isf_hdiv
26   !! * Substitutions
27#  include "do_loop_substitute.h90"
28
29CONTAINS
30
31   SUBROUTINE isf_hdiv( kt, Kmm, phdiv )
32      !!----------------------------------------------------------------------
33      !!                  ***  SUBROUTINE isf_hdiv  ***
34      !!       
35      !! ** Purpose :   update the horizontal divergence with the ice shelf contribution
36      !!                (parametrisation, explicit, ice sheet coupling conservation
37      !!                 increment)
38      !!
39      !!----------------------------------------------------------------------
40      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdiv   ! horizontal divergence
41      !!----------------------------------------------------------------------
42      INTEGER, INTENT(in) :: kt
43      INTEGER, INTENT(in) :: Kmm      !  ocean time level index
44      !
45      IF ( ln_isf ) THEN
46         !
47         ! ice shelf cavity contribution
48         IF ( ln_isfcav_mlt ) CALL isf_hdiv_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, fwfisf_cav, fwfisf_cav_b, phdiv)
49         !
50         ! ice shelf parametrisation contribution
51         IF ( ln_isfpar_mlt ) CALL isf_hdiv_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, fwfisf_par_b, phdiv)
52         !
53         ! ice sheet coupling contribution
54         IF ( ln_isfcpl .AND. kt /= 0 ) THEN
55            !
56            ! Dynamical stability at start up after change in under ice shelf cavity geometry is achieve by correcting the divergence.
57            ! This is achieved by applying a volume flux in order to keep the horizontal divergence after remapping
58            ! the same as at the end of the latest time step. So correction need to be apply at nit000 (euler time step) and
59            ! half of it at nit000+1 (leap frog time step).
60            IF ( kt == nit000   ) CALL isf_hdiv_cpl(Kmm, risfcpl_vol       , phdiv)
61            IF ( kt == nit000+1 ) CALL isf_hdiv_cpl(Kmm, risfcpl_vol*0.5_wp, phdiv)
62            !
63            ! correct divergence every time step to remove any trend due to coupling
64            ! conservation option
65            IF ( ln_isfcpl_cons ) CALL isf_hdiv_cpl(Kmm, risfcpl_cons_vol, phdiv)
66            !
67         END IF
68         !
69      END IF
70      !
71   END SUBROUTINE isf_hdiv
72
73   SUBROUTINE isf_hdiv_mlt(ktop, kbot, phtbl, pfrac, pfwf, pfwf_b, phdiv)
74      !!----------------------------------------------------------------------
75      !!                  ***  SUBROUTINE sbc_isf_div  ***
76      !!       
77      !! ** Purpose :   update the horizontal divergence with the ice shelf inflow
78      !!
79      !! ** Method  :   pfwf is positive (outflow) and expressed as kg/m2/s
80      !!                increase the divergence
81      !!
82      !! ** Action  :   phdivn   increased by the ice shelf outflow
83      !!----------------------------------------------------------------------
84      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: phdiv
85      !!----------------------------------------------------------------------
86      INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) :: ktop , kbot
87      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pfrac, phtbl
88      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pfwf , pfwf_b
89      !!----------------------------------------------------------------------
90      INTEGER  ::   ji, jj, jk   ! dummy loop indices
91      INTEGER  ::   ikt, ikb 
92      REAL(wp), DIMENSION(jpi,jpj) :: zhdiv
93      !!----------------------------------------------------------------------
94      !
95      !==   fwf distributed over several levels   ==!
96      !
97      ! compute integrated divergence correction
98      zhdiv(:,:) = 0.5_wp * ( pfwf(:,:) + pfwf_b(:,:) ) * r1_rho0 / phtbl(:,:)
99      !
100      ! update divergence at each level affected by ice shelf top boundary layer
101      DO_2D_11_11
102         ikt = ktop(ji,jj)
103         ikb = kbot(ji,jj)
104         ! level fully include in the ice shelf boundary layer
105         DO jk = ikt, ikb - 1
106            phdiv(ji,jj,jk) = phdiv(ji,jj,jk) + zhdiv(ji,jj)
107         END DO
108         ! level partially include in ice shelf boundary layer
109         phdiv(ji,jj,ikb) = phdiv(ji,jj,ikb) + zhdiv(ji,jj) * pfrac(ji,jj)
110      END_2D
111      !
112   END SUBROUTINE isf_hdiv_mlt
113
114   SUBROUTINE isf_hdiv_cpl(Kmm, pqvol, phdiv)
115      !!----------------------------------------------------------------------
116      !!                  ***  SUBROUTINE isf_hdiv_cpl  ***
117      !!       
118      !! ** Purpose :   update the horizontal divergence with the ice shelf
119      !!                coupling conservation increment
120      !!
121      !! ** Method  :   pqvol is positive (outflow) and expressed as m3/s
122      !!                increase the divergence
123      !!
124      !! ** Action  :   phdivn   increased by the ice shelf outflow
125      !!
126      !!----------------------------------------------------------------------
127      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: phdiv
128      !!----------------------------------------------------------------------
129      INTEGER,                          INTENT(in)    :: Kmm     ! ocean time level index
130      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pqvol
131      !!----------------------------------------------------------------------
132      INTEGER :: jk
133      !!----------------------------------------------------------------------
134      !
135      DO jk=1,jpk 
136         phdiv(:,:,jk) =  phdiv(:,:,jk) + pqvol(:,:,jk) * r1_e1e2t(:,:) / e3t(:,:,jk,Kmm)
137      END DO
138      !
139   END SUBROUTINE isf_hdiv_cpl
140
141END MODULE isfhdiv
Note: See TracBrowser for help on using the repository browser.