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.
agrif_dom_update.F90 in utils/tools/DOMAINcfg/src – NEMO

source: utils/tools/DOMAINcfg/src/agrif_dom_update.F90 @ 15347

Last change on this file since 15347 was 15347, checked in by jchanut, 3 years ago

#2638, set parent e3t as the maximum not the average of child grid values (without vertical remapping). Number of vertical levels was ok, though.
File size: 9.3 KB
Line 
1MODULE agrif_dom_update
2
3   USE dom_oce
4   USE agrif_parameters
5   USE agrif_profiles
6   USE agrif_recompute_scales
7 
8   IMPLICIT none
9   PRIVATE
10
11   PUBLIC agrif_update_all
12
13CONTAINS 
14
15#if defined key_agrif
16
17   SUBROUTINE agrif_update_all
18      !!----------------------------------------------------------------------
19      !!                  ***  ROUTINE agrif_update_all  ***
20      !!---------------------------------------------------------------------- 
21      !
22      INTEGER :: ind1, ind2
23
24      IF( Agrif_Root() ) return
25     
26      IF ( .NOT.ln_vert_remap ) THEN
27         CALL agrif_update_variable(bottom_level_id,procname = update_bottom_level)
28         Agrif_UseSpecialValueInUpdate = .FALSE.
29         Agrif_SpecialValueFineGrid    = 0._wp         
30         CALL agrif_update_variable(e3t_copy_id, procname = update_e3t_z) 
31         !
32      ELSE
33         Agrif_UseSpecialValueInUpdate = .FALSE.
34         Agrif_SpecialValueFineGrid    = 0._wp         
35         CALL agrif_update_variable(e3t_id, procname = update_e3t_z_cons) 
36
37         ! jc: extend update zone outside dynamical interface within sponge zone:
38         ! Use max operator this time to account for cases for which Agrif_Rho > nbghostcells
39         ind1 = CEILING(REAL(max(nbghostcells_x_w-1, nbghostcells_x_e-1), wp) / Agrif_Rhox() )
40         ind2 = CEILING(REAL(max(nbghostcells_y_s-1, nbghostcells_y_n-1), wp) / Agrif_Rhoy() )
41         CALL agrif_update_variable(e3t_copy_id, locupdate1=(/-ind1,0/), &
42                             &                   locupdate2=(/-ind2,0/),procname = update_e3t_z_cons)
43      ENDIF
44      Agrif_UseSpecialValueInUpdate = .FALSE.
45      !
46      ! Update vertical scale factors at U, V and F-points:
47      CALL Agrif_ChildGrid_To_ParentGrid()
48      CALL agrif_recompute_scalefactors
49      CALL Agrif_ParentGrid_To_ChildGrid()
50      !   
51   END SUBROUTINE agrif_update_all
52
53   SUBROUTINE update_bottom_level( ptab, i1, i2, j1, j2, before)
54      !!----------------------------------------------------------------------
55      !!       ***  ROUTINE update_bottom_level  ***
56      !!---------------------------------------------------------------------- 
57      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
58      REAL, DIMENSION(i1:i2,j1:j2)    , INTENT(inout) ::   ptab
59      LOGICAL                         , INTENT(in   ) ::   before
60      !
61      !!----------------------------------------------------------------------
62      !
63      IF( before) THEN
64         ptab(i1:i2,j1:j2) = mbkt(i1:i2,j1:j2)*ssmask(i1:i2,j1:j2)
65      ELSE
66         mbkt(i1:i2,j1:j2) = nint(ptab(i1:i2,j1:j2))
67         
68         WHERE ( mbkt(i1:i2,j1:j2) .EQ. 0 )
69            ssmask(i1:i2,j1:j2) = 0._wp
70            mbkt(i1:i2,j1:j2)   = 1 
71         ELSEWHERE
72            ssmask(i1:i2,j1:j2) = 1._wp
73         END WHERE
74      ENDIF
75      !
76   END SUBROUTINE update_bottom_level
77
78   SUBROUTINE update_e3t_z( tabres, i1, i2, j1, j2, k1, k2, before )
79      !!---------------------------------------------
80      !!           *** update_e3t_z ***
81      !!---------------------------------------------
82      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
83      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
84      LOGICAL, INTENT(in) :: before
85      !!
86      INTEGER :: ji, jj, jk
87      !!---------------------------------------------
88      !
89      IF (before) THEN
90         DO jk=k1,k2
91            DO jj=j1,j2
92               DO ji=i1,i2
93                  IF ( (ssmask(ji,jj) /=0._wp).AND.(mbkt(ji,jj).GE.jk) ) THEN
94                     tabres(ji,jj,jk) = e3t_0(ji,jj,jk)
95                  ELSE
96                     tabres(ji,jj,jk) = 0._wp
97                  ENDIF
98               END DO
99            END DO
100         END DO
101      ELSE
102         DO jk=1,jpk
103            DO jj=j1,j2
104               DO ji=i1,i2
105                  IF ( ( mbkt(ji,jj).GE.jk ).AND.(ssmask(ji,jj)==1._wp) ) THEN
106                     e3t_0(ji,jj,jk) = MAX(tabres(ji,jj,jk),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat))
107                 !    e3t_0(ji,jj,jk) = tabres(ji,jj,jk)
108                  ELSE
109                     e3t_0(ji,jj,jk) = e3t_1d(jk)
110                  ENDIF
111               END DO
112            END DO
113         END DO
114         !
115      ENDIF
116      !
117   END SUBROUTINE update_e3t_z
118
119   SUBROUTINE update_e3t_z_cons( tabres, i1, i2, j1, j2, k1, k2, before )
120      !!---------------------------------------------
121      !!           *** update_e3t_z_cons ***
122      !!---------------------------------------------
123      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
124      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
125      LOGICAL, INTENT(in) :: before
126      !!
127      INTEGER :: ji, jj, jk, ik
128      REAL(wp) :: zhmin, zdepth, zdepwp, ze3tp, ze3wp
129      !!---------------------------------------------
130      !
131      IF (before) THEN
132         DO jk = k1, k2-1
133            DO jj = j1, j2
134               DO ji = i1, i2
135                   IF ( (ssmask(ji,jj) /=0._wp).AND.( mbkt(ji,jj) .GE. jk ) ) THEN
136                      tabres(ji,jj,jk) = e3t_0(ji,jj,jk)
137                   ELSE
138                      tabres(ji,jj,jk) = 0._wp
139                   endif
140               END DO
141            END DO
142         END DO
143         tabres(i1:i2,j1:j2,k2) = ssmask(i1:i2,j1:j2) ! To get fractional area
144      ELSE
145         IF( rn_hmin < 0._wp ) THEN   
146            ik = - INT( rn_hmin )
147         ELSE                         
148            ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )
149         ENDIF
150         zhmin = gdepw_1d(ik+1)
151
152         ! Compute child bathymetry:
153         bathy(i1:i2,j1:j2) = 0._wp
154         DO jk=k1,k2-1   
155            bathy(i1:i2,j1:j2) = bathy(i1:i2,j1:j2) + tabres(i1:i2,j1:j2,jk)
156         END DO
157         WHERE( bathy(i1:i2,j1:j2) == 0._wp )   ;   mbathy(i1:i2,j1:j2) = 0       
158         ELSE WHERE                             ;   mbathy(i1:i2,j1:j2) = jpkm1 
159         END WHERE
160
161         DO jk = jpkm1, 1, -1
162            zdepth = gdepw_1d(jk) ! + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
163            WHERE( 0._wp < bathy(i1:i2,j1:j2) .AND. bathy(i1:i2,j1:j2) <= zdepth )   mbathy(i1:i2,j1:j2) = jk-1
164         END DO
165
166         ! Scale factors and depth at T- and W-points
167         DO jk = 1, jpk 
168            gdept_0(i1:i2,j1:j2,jk) = gdept_1d(jk)
169            gdepw_0(i1:i2,j1:j2,jk) = gdepw_1d(jk)
170            e3t_0  (i1:i2,j1:j2,jk) = e3t_1d  (jk)
171            e3w_0  (i1:i2,j1:j2,jk) = e3w_1d  (jk)
172         END DO
173         ! Scale factors and depth at T- and W-points
174         DO jj = j1, j2
175            DO ji = i1, i2 
176               ik = mbathy(ji,jj)
177               IF( ik > 0 ) THEN               ! ocean point only
178                  ! max ocean level case
179                  IF( ik == jpkm1 ) THEN
180                     zdepwp = bathy(ji,jj)
181                     ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
182                     ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
183                     e3t_0(ji,jj,ik  ) = ze3tp
184                     e3t_0(ji,jj,ik+1) = ze3tp
185                     e3w_0(ji,jj,ik  ) = ze3wp
186                     e3w_0(ji,jj,ik+1) = ze3tp
187                     gdepw_0(ji,jj,ik+1) = zdepwp
188                     gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
189                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
190                     !
191                  ELSE                         ! standard case
192                     IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN
193                        gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
194                     ELSE                                       
195                        gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
196                     ENDIF
197                     gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )           &
198                        &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )           &
199                        &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
200                     e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik))           &
201                        &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik)) 
202                     e3w_0(ji,jj,ik) =                                                                   & 
203                        &      0.5_wp * (gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
204                        &             * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
205                     !       ... on ik+1
206                     e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
207                     e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
208                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
209                  ENDIF
210               ENDIF
211            END DO
212         END DO
213         !
214         DO jj=j1,j2
215            DO ji=i1,i2
216               bathy(ji,jj) = SUM(e3t_0(ji,jj,1:mbkt(ji,jj) ) ) 
217            END DO
218         END DO
219         !
220         WHERE ( ( mbathy(i1:i2,j1:j2) .EQ. 0 )  & 
221           & .OR.(tabres(i1:i2,j1:j2,k2)<0.5_wp) &
222           & .OR.(bathy(i1:i2,j1:j2)<zhmin) )
223            ssmask(i1:i2,j1:j2) = 0._wp
224            mbathy(i1:i2,j1:j2) = 0 
225         ELSEWHERE
226            ssmask(i1:i2,j1:j2) = 1._wp
227         END WHERE
228         mbkt(i1:i2,j1:j2) = MAX( mbathy(i1:i2,j1:j2), 1 )
229      ENDIF
230      !
231   END SUBROUTINE update_e3t_z_cons
232     
233#else
234   SUBROUTINE agrif_update_all
235   END SUBROUTINE agrif_update_all
236#endif
237
238END MODULE agrif_dom_update
Note: See TracBrowser for help on using the repository browser.