source: branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscpldiv.F90 @ 9813

Last change on this file since 9813 was 9813, checked in by antsia, 2 years ago

delete iscplhsb, add iscpldiv and make the code readable

File size: 8.2 KB
Line 
1MODULE iscpldiv
2   !!======================================================================
3   !!                       ***  MODULE  iscpldiv***
4   !! Divergence correction : compute and save the divergence correction due to ice sheet/ocean coupling
5   !!=====================================================================
6   !! History :  NEMO  ! 2018-06 A. Siahaan: original
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   iscpl_div_corr     : compute the divergence correction
11   !!   iscpl_div          : update hdivn with the divergence correction
12   !!   iscpl_div_tra      : update tsa with contribution from the divergence correction
13   !!----------------------------------------------------------------------
14   USE dom_oce         ! ocean space and time domain
15   USE oce             ! global tra/dyn variable
16   USE in_out_manager  ! I/O manager
17   USE lib_mpp         ! MPP library
18   USE lib_fortran     ! MPP library
19   USE wrk_nemo        ! Memory allocation
20   USE lbclnk          !
21   USE iscplini
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC iscpl_div
27   PUBLIC iscpl_div_tra
28   PUBLIC iscpl_div_corr
29   !!
30   !! * Substitutions 
31#  include "domzgr_substitute.h90"
32#  include "vectopt_loop_substitute.h90"
33   !!----------------------------------------------------------------------
34   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
35   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
36   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
37   !!----------------------------------------------------------------------
38CONTAINS
39
40   SUBROUTINE iscpl_div( phdivn )
41      !!----------------------------------------------------------------------
42      !!                  ***  ROUTINE iscpl_div  ***
43      !!
44      !! ** Purpose :   update the horizontal divergence
45      !!
46      !! ** Method  :
47      !!                CAUTION : iscpl is positive (inflow) and expressed in
48      !m/s
49      !!
50      !! ** Action  :   phdivn   increase by the iscpl correction term
51      !!----------------------------------------------------------------------
52      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
53      !!
54      !!----------------------------------------------------------------------
55      !
56      IF(lwp) WRITE(numout,*) 'adding the divergence correction contribution at nit000'
57
58      phdivn(:,:,:) =  phdivn(:,:,:) + rhdivdiff(:,:,:)
59
60   END SUBROUTINE iscpl_div
61
62   SUBROUTINE iscpl_div_tra( ptsa )
63      !!----------------------------------------------------------------------
64      !!                  ***  ROUTINE iscpl_div_tra  ***
65      !!
66      !! ** Purpose :   update the active tracer tendency
67      !!
68      !! ** Method  :
69      !!                CAUTION : iscpl is positive (inflow) and expressed in
70      !m/s
71      !!
72      !! ** Action  :   ptsa   increase by the iscpl correction term
73      !!----------------------------------------------------------------------
74      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   ptsa   ! active tracer tendency
75      !!
76      !!----------------------------------------------------------------------
77      !
78      IF (lwp) WRITE(numout,*) 'add tendency tsa due to re-meshing at kt = nit000'
79
80      ptsa(:,:,:,1) = ptsa(:,:,:,1) + rhdivdiff_trac(:,:,:,1)*tmask(:,:,:)
81      ptsa(:,:,:,2) = ptsa(:,:,:,2) + rhdivdiff_trac(:,:,:,2)*tmask(:,:,:)
82
83   END SUBROUTINE iscpl_div_tra
84
85
86
87   SUBROUTINE iscpl_div_corr(ptmask_b, pe3t_b, pe3u_b, pe3v_b)
88      !!----------------------------------------------------------------------
89      !!                  ***  ROUTINE iscpl_div_corr  ***
90      !!
91      !! ** Purpose :   compute the divergence correction
92      !!
93      !! ** Method  :
94      !!               compute the difference of hdivn*e3t_n between before and after remeshing   
95      !!               transfer the missing hdvn*e3t_n from top of isf to the new mikt
96      !!               transfer the missing hdvn*e3t_n from bottom of isf to the new mbkt
97      !!----------------------------------------------------------------------
98      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b !! mask before
99      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
100      !!
101      INTEGER :: ji, jj, jk, jtop         !! loop index
102      !!
103      REAL(wp):: zvol, zuflux_sum, zvflux_sum
104
105      REAL(wp), DIMENSION(:,:,:  ), POINTER :: zhdiv, zhdiv2
106      !!----------------------------------------------------------------------
107
108      !! allocate variables
109      CALL wrk_alloc(jpi,jpj,jpk, zhdiv, zhdiv2)
110
111      !-----initialise the divergence and tracer correction------
112      zhdiv(:,:,:) = 0._wp
113      rhdivdiff(:,:,:) = 0._wp
114      rhdivdiff_trac(:,:,:,:) = 0._wp     
115      zhdiv2(:,:,:) = 0._wp
116
117
118     !calculate the hdivn*e3tn with the last values before re-meshing (should only be different from the real hdivn for cavity only)
119     !hdivb and rot should remain the same as the previous timestep bcause it will be used in dyn_ldf, so we did not use div_cur to calculate rhdivdiff as we should only use div_cur once for div_cur(nit000) in order not to update hdivb and rot
120
121!----------------------------------------------hdivn before remeshing-------------------------------------
122   DO jj = 2, jpjm1
123         DO ji = 2, jpim1
124            DO jk = 1, jpk                                 ! Horizontal slab
125               zhdiv2(ji,jj,jk) =   &
126                  (  e2u(ji,jj)*pe3u_b(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  )*pe3u_b(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)       &
127                   + e1v(ji,jj)*pe3v_b(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1)*pe3v_b(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  )    &
128                  / ( e1t(ji,jj) * e2t(ji,jj) )*ptmask_b(ji,jj,jk)
129            END DO
130         END DO
131      ENDDO
132
133! velocity with the new mask
134      un(:,:,:)=un(:,:,:)*umask(:,:,:)
135      vn(:,:,:)=vn(:,:,:)*vmask(:,:,:)
136
137
138      !calling divcur to update hdivn (ingat: hdivb = hdivn(before re-meshing), yet that hdivn was not using the recent values, as it was not updated at the end of timestep  )
139      !CALL div_cur(nit000)
140!----------------------------------------------hdivn after remeshing-------------------------------------
141      DO jj = 2, jpjm1
142         DO ji = 2, jpim1
143            DO jk = 1, jpk                                 ! Horizontal slab
144                 zhdiv(ji,jj,jk) =   &
145                  (  e2u(ji,jj)*fse3u_n(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  )*fse3u_n(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)       &
146                   + e1v(ji,jj)*fse3v_n(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1)*fse3v_n(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  )    &
147                  / ( e1t(ji,jj) * e2t(ji,jj) )*tmask(ji,jj,jk)
148            END DO
149         END DO
150      ENDDO
151
152      zhdiv2(:,:,:) = zhdiv2(:,:,:) - zhdiv(:,:,:)  ! non-cavity will be = 0
153      rhdivdiff(:,:,:) =  zhdiv2(:,:,:) 
154
155      !----transfers the divergence of closed cells to the new top cell------
156      !until this point, rhdivdiff has the same unit as hdivn*e3t
157      DO jj = 2, jpjm1
158        DO ji = 2, jpim1
159          jk = mikt(ji,jj) 
160         IF (jk > 1 .AND. ssmask(ji,jj) == 1 ) THEN
161            rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + SUM(rhdivdiff(ji,jj,1:jk-1)*ptmask_b(ji,jj,1:jk-1)) 
162         ENDIF
163
164         jk = mbkt(ji,jj)
165         IF (ssmask(ji,jj) == 1 ) THEN       
166            rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + SUM(rhdivdiff(ji,jj,jk+1:jpk)*ptmask_b(ji,jj,jk+1:jpk)) 
167        ENDIF
168
169        END DO
170      ENDDO
171 
172
173      !----cancel all the divergence of closed cells with tmask, then convert it to hdivn by dividing with e3t_n (it is now safe to do this)
174      ! also cancelling totally closed cavity column
175      rhdivdiff(:,:,:) = rhdivdiff(:,:,:)*tmask(:,:,:)/fse3t_n(:,:,:)
176      !tracer correction trend associated with the divergence correction; opposite sign of divergence
177      rhdivdiff_trac(:,:,:,1) = - rhdivdiff(:,:,:) * tsb(:,:,:,1)*tmask(:,:,:)
178      rhdivdiff_trac(:,:,:,2) = - rhdivdiff(:,:,:) * tsb(:,:,:,2)*tmask(:,:,:)
179
180      CALL lbc_lnk(rhdivdiff,   'T',1._wp)
181      CALL lbc_lnk(rhdivdiff_trac(:,:,:,1),   'T',1._wp)
182      CALL lbc_lnk(rhdivdiff_trac(:,:,:,2),   'T',1._wp)
183
184     CALL wrk_dealloc(jpi,jpj,jpk, zhdiv, zhdiv2 )
185
186   END SUBROUTINE iscpl_div_corr
187
188END MODULE iscpldiv
Note: See TracBrowser for help on using the repository browser.