- Timestamp:
- 2016-10-29T01:21:05+02:00 (8 years ago)
- Location:
- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/NST_SRC
- Files:
-
- 3 edited
-
agrif_lim3_interp.F90 (modified) (5 diffs)
-
agrif_lim3_update.F90 (modified) (2 diffs)
-
agrif_user.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90
r7069 r7158 24 24 USE ice 25 25 USE agrif_ice 26 26 27 27 IMPLICIT NONE 28 28 PRIVATE … … 38 38 CONTAINS 39 39 40 SUBROUTINE agrif_interp_lim3( cd_type )40 SUBROUTINE agrif_interp_lim3( cd_type, kiter, kitermax ) 41 41 !!----------------------------------------------------------------------- 42 42 !! *** ROUTINE agrif_rhg_lim3 *** 43 43 !! 44 44 !! ** Method : simple call to atomic routines using stored values to 45 !! fill the boundaries depending of the position of the point and 46 !! computing factor for time interpolation 47 !!----------------------------------------------------------------------- 48 CHARACTER(len=1), INTENT( in ) :: cd_type 49 !! 45 !! fill the boundaries depending of the position of the point and 46 !! computing factor for time interpolation 47 !!----------------------------------------------------------------------- 48 CHARACTER(len=1), INTENT( in ) :: cd_type 49 INTEGER , INTENT( in ), OPTIONAL :: kiter, kitermax 50 !! 50 51 REAL(wp) :: zbeta 51 52 !!----------------------------------------------------------------------- … … 53 54 IF( Agrif_Root() ) RETURN 54 55 ! 55 zbeta = REAL(lim_nbstep) / ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) ) 56 ! 57 ! clem: calledweight = zbeta(1/3;2/3;1) => 2/3*var1+1/3*var2 puis 1/3;2/3 puis 0;1 56 IF( PRESENT( kiter ) ) THEN ! interpolation at the child sub-time step (for ice rheology) 57 zbeta = ( REAL(lim_nbstep) - REAL(kitermax - kiter) / REAL(kitermax) ) / & 58 & ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) ) 59 ELSE ! interpolation at the child time step 60 zbeta = REAL(lim_nbstep) / ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) ) 61 ENDIF 62 ! 58 63 Agrif_SpecialValue=-9999. 59 64 Agrif_UseSpecialValue = .TRUE. … … 126 131 127 132 128 SUBROUTINE interp_tra_ice( ptab, i1, i2, j1, j2, k1, k2, before )133 SUBROUTINE interp_tra_ice( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir ) 129 134 !!----------------------------------------------------------------------- 130 135 !! *** ROUTINE interp_tra_ice *** … … 137 142 INTEGER , INTENT(in) :: i1, i2, j1, j2, k1, k2 138 143 LOGICAL , INTENT(in) :: before 139 !! 140 INTEGER :: jk, jl, jm 144 INTEGER , INTENT(in) :: nb, ndir 145 !! 146 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztab 147 INTEGER :: ji, jj, jk, jl, jm 148 INTEGER :: imin, imax, jmin, jmax 149 REAL(wp) :: zrhox, z1, z2, z3, z4, z5, z6, z7 150 LOGICAL :: western_side, eastern_side, northern_side, southern_side 151 141 152 !!----------------------------------------------------------------------- 142 153 ! clem: pkoi on n'utilise pas les quantités intégrées ici => before: * e12t ; after: * r1_e12t / rhox / rhoy 143 154 ! a priori c'est ok comme ca (cf ce qui est fait dans l'ocean). Je ne sais pas pkoi ceci dit 144 155 ALLOCATE( ztab(SIZE(a_i_b,1),SIZE(a_i_b,2),SIZE(ptab,3)) ) 156 145 157 IF( before ) THEN ! parent grid 146 158 jm = 1 147 159 DO jl = 1, jpl 148 ptab( :,:,jm) = a_i_b (i1:i2,j1:j2,jl) ; jm = jm + 1149 ptab( :,:,jm) = v_i_b (i1:i2,j1:j2,jl) ; jm = jm + 1150 ptab( :,:,jm) = v_s_b (i1:i2,j1:j2,jl) ; jm = jm + 1151 ptab( :,:,jm) = smv_i_b(i1:i2,j1:j2,jl) ; jm = jm + 1152 ptab( :,:,jm) = oa_i_b (i1:i2,j1:j2,jl) ; jm = jm + 1160 ptab(i1:i2,j1:j2,jm) = a_i_b (i1:i2,j1:j2,jl) ; jm = jm + 1 161 ptab(i1:i2,j1:j2,jm) = v_i_b (i1:i2,j1:j2,jl) ; jm = jm + 1 162 ptab(i1:i2,j1:j2,jm) = v_s_b (i1:i2,j1:j2,jl) ; jm = jm + 1 163 ptab(i1:i2,j1:j2,jm) = smv_i_b(i1:i2,j1:j2,jl) ; jm = jm + 1 164 ptab(i1:i2,j1:j2,jm) = oa_i_b (i1:i2,j1:j2,jl) ; jm = jm + 1 153 165 DO jk = 1, nlay_s 154 ptab( :,:,jm) = e_s_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1166 ptab(i1:i2,j1:j2,jm) = e_s_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1 155 167 ENDDO 156 168 DO jk = 1, nlay_i 157 ptab( :,:,jm) = e_i_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1169 ptab(i1:i2,j1:j2,jm) = e_i_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1 158 170 ENDDO 159 171 ENDDO 160 !!ptab(:,:,jm) = ato_i(i1:i2,j1:j2)161 172 162 173 DO jk = k1, k2 163 WHERE( tmask(i1:i2,j1:j2,1) == 0. ) ptab( :,:,jk) = -9999.174 WHERE( tmask(i1:i2,j1:j2,1) == 0. ) ptab(i1:i2,j1:j2,jk) = -9999. 164 175 ENDDO 165 176 166 177 ELSE ! child grid 178 !! ==> The easiest interpolation is the following commented lines 179 !! jm = 1 180 !! DO jl = 1, jpl 181 !! a_i (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 182 !! v_i (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 183 !! v_s (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 184 !! smv_i(i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 185 !! oa_i (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 186 !! DO jk = 1, nlay_s 187 !! e_s(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 188 !! ENDDO 189 !! DO jk = 1, nlay_i 190 !! e_i(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 191 !! ENDDO 192 !! ENDDO 193 194 !! ==> this is a more complex interpolation since we mix solutions over a couple of grid points 195 !! it is advised to use it for fields modified by high order schemes (e.g. advection UM5...) 196 ! record ztab 167 197 jm = 1 168 198 DO jl = 1, jpl 169 a_i (i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1170 v_i (i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1171 v_s (i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1172 smv_i(i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1173 oa_i (i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1199 ztab(:,:,jm) = a_i_b (:,:,jl) ; jm = jm + 1 200 ztab(:,:,jm) = v_i_b (:,:,jl) ; jm = jm + 1 201 ztab(:,:,jm) = v_s_b (:,:,jl) ; jm = jm + 1 202 ztab(:,:,jm) = smv_i_b(:,:,jl) ; jm = jm + 1 203 ztab(:,:,jm) = oa_i_b (:,:,jl) ; jm = jm + 1 174 204 DO jk = 1, nlay_s 175 e_s(i1:i2,j1:j2,jk,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1205 ztab(:,:,jm) = e_s_b(:,:,jk,jl) ; jm = jm + 1 176 206 ENDDO 177 207 DO jk = 1, nlay_i 178 e_i(i1:i2,j1:j2,jk,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1208 ztab(:,:,jm) = e_i_b(:,:,jk,jl) ; jm = jm + 1 179 209 ENDDO 180 210 ENDDO 181 !!ato_i(i1:i2,j1:j2) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) 211 ! 212 ! borders of the domain 213 western_side = (nb == 1).AND.(ndir == 1) ; eastern_side = (nb == 1).AND.(ndir == 2) 214 southern_side = (nb == 2).AND.(ndir == 1) ; northern_side = (nb == 2).AND.(ndir == 2) 215 ! 216 ! spatial smoothing 217 zrhox = Agrif_Rhox() 218 z1 = ( zrhox - 1. ) * 0.5 219 z3 = ( zrhox - 1. ) / ( zrhox + 1. ) 220 z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. ) 221 z7 = - ( zrhox - 1. ) / ( zrhox + 3. ) 222 z2 = 1. - z1 223 z4 = 1. - z3 224 z5 = 1. - z6 - z7 225 ! 226 ! Remove corners 227 imin = i1 ; imax = i2 ; jmin = j1 ; jmax = j2 228 IF( (nbondj == -1) .OR. (nbondj == 2) ) jmin = 3 229 IF( (nbondj == +1) .OR. (nbondj == 2) ) jmax = nlcj-2 230 IF( (nbondi == -1) .OR. (nbondi == 2) ) imin = 3 231 IF( (nbondi == +1) .OR. (nbondi == 2) ) imax = nlci-2 232 233 ! smoothed fields 234 IF( eastern_side ) THEN 235 ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:) 236 DO jj = jmin, jmax 237 rswitch = 0. 238 IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1. 239 ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:) & 240 & + umask(nlci-2,jj,1) * & 241 & ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:) + z3 * ztab(nlci-2,jj,:) ) & 242 & + rswitch * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) ) 243 ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1) 244 END DO 245 ENDIF 246 ! 247 IF( northern_side ) THEN 248 ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:) 249 DO ji = imin, imax 250 rswitch = 0. 251 IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1. 252 ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:) & 253 & + vmask(ji,nlcj-2,1) * & 254 & ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:) + z3 * ztab(ji,nlcj-2,:) ) & 255 & + rswitch * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) ) 256 ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1) 257 END DO 258 END IF 259 ! 260 IF( western_side) THEN 261 ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:) 262 DO jj = jmin, jmax 263 rswitch = 0. 264 IF( u_ice(2,jj) > 0._wp ) rswitch = 1. 265 ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:) & 266 & + umask(2,jj,1) * & 267 & ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) & 268 & + rswitch * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) ) 269 ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1) 270 END DO 271 ENDIF 272 ! 273 IF( southern_side ) THEN 274 ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:) 275 DO ji = imin, imax 276 rswitch = 0. 277 IF( v_ice(ji,2) > 0._wp ) rswitch = 1. 278 ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:) & 279 & + vmask(ji,2,1) * & 280 & ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) & 281 & + rswitch * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) ) 282 ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1) 283 END DO 284 END IF 285 ! 286 ! Treatment of corners 287 IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) ) ztab(nlci-1,2,:) = ptab(nlci-1,2,:) ! East south 288 IF( (eastern_side) .AND. ((nbondj == 1).OR.(nbondj == 2)) ) ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north 289 IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) ) ztab(2,2,:) = ptab(2,2,:) ! West south 290 IF( (western_side) .AND. ((nbondj == 1).OR.(nbondj == 2)) ) ztab(2,nlcj-1,:) = ptab(2,nlcj-1,:) ! West north 291 292 ! retrieve ice tracers 293 jm = 1 294 DO jl = 1, jpl 295 a_i (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 296 v_i (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 297 v_s (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 298 smv_i(i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 299 oa_i (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 300 DO jk = 1, nlay_s 301 e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 302 ENDDO 303 DO jk = 1, nlay_i 304 e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1 305 ENDDO 306 ENDDO 182 307 183 308 ENDIF 309 310 DEALLOCATE( ztab ) 184 311 ! 185 312 END SUBROUTINE interp_tra_ice -
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/NST_SRC/agrif_lim3_update.F90
r6746 r7158 52 52 !!---------------------------------------------------------------------- 53 53 ! 54 IF( ( Agrif_NbStepint() .NE. (Agrif_irhot()-1) ) .AND. (kt /= 0) ) RETURN ! do not update if nb of child time steps differ from time refinement 55 ! i.e. update only at the parent time step 54 ! IF( ( MOD( kt-nit000, Agrif_irhot() * Agrif_Parent(nn_fsbc) ) /=0 ) .AND. (kt /= 0) ) THEN 55 ! PRINT *, 'clem NOT udpate, kt=',kt,Agrif_NbStepint() 56 ! ELSE 57 ! PRINT *, 'clem UPDATE, kt=',kt,Agrif_NbStepint() 58 ! ENDIF 56 59 60 !! clem: I think the update should take place each time the ocean sees the surface forcings (but maybe I am wrong and we should update every rhot time steps) 61 IF( ( MOD( kt-nit000, Agrif_irhot() * Agrif_Parent(nn_fsbc) ) /=0 ) .AND. (kt /= 0) ) RETURN ! do not update if nb of child time steps differ from time refinement 62 ! i.e. update only at the parent time step 63 !! clem: this condition is clearly wrong if nn_fsbc/=1 (==> Agrif_NbStepint /= (Agrif_irhot()-1) all the time) 64 !!IF( ( Agrif_NbStepint() .NE. (Agrif_irhot()-1) ) .AND. (kt /= 0) ) RETURN 65 57 66 Agrif_UseSpecialValueInUpdate = .TRUE. 58 67 Agrif_SpecialValueFineGrid = -9999. … … 60 69 IF( MOD(nbcline,nbclineupdate) == 0) THEN ! update the whole basin at each nbclineupdate (=nn_cln_update) baroclinic parent time steps 61 70 ! nbcline is incremented (+1) at the end of each parent time step from 0 (1st time step) 62 ! clem: j'ai l'impression qu'il y a un decalage de 1 mais selon rachid c ok63 71 CALL Agrif_Update_Variable( tra_ice_id , procname = update_tra_ice ) 64 72 CALL Agrif_Update_Variable( u_ice_id , procname = update_u_ice ) 65 73 CALL Agrif_Update_Variable( v_ice_id , procname = update_v_ice ) 66 ELSE ! update only the boundaries 67 ! defined par locupdate 74 ELSE ! update only the boundaries defined par locupdate 68 75 CALL Agrif_Update_Variable( tra_ice_id , locupdate=(/0,2/), procname = update_tra_ice ) 69 76 CALL Agrif_Update_Variable( u_ice_id , locupdate=(/0,1/), procname = update_u_ice ) -
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
r7060 r7158 241 241 WRITE(cl_check2,*) NINT(rdt) 242 242 WRITE(cl_check3,*) NINT(Agrif_Parent(rdt)/Agrif_Rhot()) 243 CALL ctl_ warn( 'incompatible time step between grids', &243 CALL ctl_stop( 'incompatible time step between ocean grids', & 244 244 & 'parent grid value : '//cl_check1 , & 245 245 & 'child grid value : '//cl_check2 , & 246 & 'value on child grid will be changed to : '//cl_check3 ) 247 rdt=Agrif_Parent(rdt)/Agrif_Rhot() 246 & 'value on child grid should be changed to : '//cl_check3 ) 248 247 ENDIF 249 248 … … 565 564 !!---------------------------------------------------------------------- 566 565 USE Agrif_Util 566 USE sbc_oce, ONLY : nn_fsbc ! clem: necessary otherwise Agrif_Parent(nn_fsbc) = nn_fsbc 567 567 USE ice 568 568 USE agrif_ice … … 579 579 CALL agrif_declare_var_lim3 580 580 581 ! clem: reset nn_fsbc(child) to rhot if rhot * nn_fsbc(parent) /= N * nn_fsbc(child) with N being integer 581 ! Controls (clem) 582 ! stop if rhot * nn_fsbc(parent) /= N * nn_fsbc(child) with N being integer 582 583 IF( MOD( Agrif_irhot() * Agrif_Parent(nn_fsbc), nn_fsbc ) /= 0 ) THEN 583 nn_fsbc = Agrif_irhot() 584 CALL ctl_warn ('rhot * nn_fsbc(parent) /= N * nn_fsbc(child), therefore nn_fsbc(child) is set to rhot') 585 WRITE(*,*) 'New nn_fsbc(child) = ', nn_fsbc 584 CALL ctl_stop('rhot * nn_fsbc(parent) /= N * nn_fsbc(child), therefore nn_fsbc(child) should be set to 1 or nn_fsbc(parent)') 586 585 ENDIF 587 586 588 ! clem: reset update frequency if different from nn_fsbc 589 IF( nbclineupdate /= nn_fsbc ) THEN 590 nbclineupdate = nn_fsbc 591 CALL ctl_warn ('With ice model on child grid, nc_cln_update is set to nn_fsbc') 592 ENDIF 587 ! stop if update frequency is different from nn_fsbc 588 IF( nbclineupdate > nn_fsbc ) CALL ctl_stop('With ice model on child grid, nn_cln_update should be set to 1 or nn_fsbc') 589 593 590 594 591 ! First Interpolations (using "after" ice subtime step => lim_nbstep=1) … … 603 600 !---------------------- 604 601 CALL agrif_update_lim3(0) 602 605 603 ! 606 604 END SUBROUTINE Agrif_InitValues_cont_lim3 … … 613 611 !!---------------------------------------------------------------------- 614 612 USE Agrif_Util 615 USE Agrif_ice !clem useless ?616 613 USE ice 617 614 … … 703 700 WRITE(cl_check2,*) rdt 704 701 WRITE(cl_check3,*) rdt*Agrif_Rhot() 705 CALL ctl_ warn( 'incompatible time step between grids', &702 CALL ctl_stop( 'incompatible time step between grids', & 706 703 & 'parent grid value : '//cl_check1 , & 707 704 & 'child grid value : '//cl_check2 , & 708 & 'value on child grid willbe changed to &705 & 'value on child grid should be changed to & 709 706 & :'//cl_check3 ) 710 rdt=rdt*Agrif_Rhot()711 707 ENDIF 712 708
Note: See TracChangeset
for help on using the changeset viewer.
