- Timestamp:
- 2014-01-28T14:55:35+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90
r4356 r4375 49 49 INTEGER :: ji, jj, jk, jk1 ! dummy loop indices 50 50 INTEGER :: zflag, z2dt ! local scalar 51 REAL(wp) :: zflxp, zflxn ! local scalars52 51 REAL(wp) :: zcoef, zdep1, zdep2 ! local scalars 53 52 REAL(wp) :: ztmp ! local scalars 53 REAL(wp), POINTER, DIMENSION(:,:) :: zflxp, zflxn ! specific 2D workspace 54 54 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu, zflxv ! specific 2D workspace 55 55 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu1, zflxv1 ! specific 2D workspace 56 REAL(wp), POINTER, DIMENSION(:,:) :: uwdlmt, vwdlmt !: W/D flux limiters 56 REAL(wp), POINTER, DIMENSION(:,:) :: wdlmt !: W/D flux limiters 57 57 58 !!---------------------------------------------------------------------- 58 59 ! … … 62 63 IF(ln_wad) THEN 63 64 64 CALL wrk_alloc( jpi, jpj, zflx u, zflxv, zflxu1, zflxv1, uwdlmt, vwdlmt)65 CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1, wdlmt ) 65 66 ! 66 67 … … 71 72 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 72 73 73 74 zflxp(:,:) = 0._wp 75 zflxn(:,:) = 0._wp 74 76 zflxu(:,:) = 0._wp 75 77 zflxv(:,:) = 0._wp 76 uwdlmt(:,:) = 1._wp 77 vwdlmt(:,:) = 1._wp 78 78 wdlmt(:,:) = 1._wp 79 79 80 80 ! Horizontal Flux in u and v direction … … 91 91 zflxv(:,:) = zflxv(:,:) * e1v(:,:) 92 92 93 DO jk1 = 1, 1094 95 zflag = 0 ! flag marking that if we need further iteration of W/D limiters?96 97 zflxu1(:,:) = zflxu(:,:) * uwdlmt(:,:)98 zflxv1(:,:) = zflxv(:,:) * vwdlmt(:,:)99 100 101 93 DO jj = 2, jpjm1 102 94 DO ji = fs_2, fs_jpim1 ! vector opt. 103 zflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) + & 104 & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,jj-1), 0._wp) 105 zflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) + & 106 & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,jj-1), 0._wp) 107 95 108 96 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 109 110 !IF( .NOT. AGRIF_Root() ) THEN !Wetting and Drying but not for AGRIF 111 IF (((nbondi == 1).OR.(nbondi == 2)).AND.(ji == nlci-1)) CYCLE ! east 112 IF (((nbondi == -1).OR.(nbondi == 2)).AND.(ji == 2 )) CYCLE ! west 113 IF (((nbondj == 1).OR.(nbondj == 2)).AND.(jj == nlcj-1)) CYCLE ! north 114 IF (((nbondj == -1).OR.(nbondj == 2)).AND.(jj == 2 )) CYCLE ! south 115 !END IF 116 117 118 ztmp = e1t(ji,jj) * e2t(ji,jj) 119 120 zdep1 = (zflxp + zflxn) * z2dt / ztmp 97 98 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & 99 & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji, jj-1), 0._wp) 100 zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj), 0._wp) + & 101 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 102 121 103 zdep2 = bathy(ji,jj) + sshb(ji,jj) - rn_wadmin 122 123 104 IF(zdep2 < 0._wp) THEN 124 !WRITE(numout,*) 'depth less than minimum depth, cell(ji,jj):', ji,jj125 105 zdep2 = 0._wp 126 106 sshb(ji,jj) = rn_wadmin - bathy(ji,jj) 127 107 END IF 108 ENDDO 109 END DO 110 111 112 !! start limiter iteration 113 DO jk1 = 1, nn_waditr 128 114 129 IF(zdep1 > zdep2) THEN 130 zflag = 1 131 zcoef = ( ( zdep2 + rn_wadmine ) * ztmp + zflxn * z2dt ) / ( zflxp * z2dt ) 132 IF(zflxu1(ji, jj) >= 0._wp) uwdlmt(ji,jj ) = zcoef 133 IF(zflxu1(ji-1,jj) < 0._wp) uwdlmt(ji-1,jj) = zcoef 134 IF(zflxv1(ji, jj) >= 0._wp) vwdlmt(ji,jj ) = zcoef 135 IF(zflxv1(ji,jj-1) < 0._wp) vwdlmt(ji,jj-1) = zcoef 136 END IF 137 END DO ! ji loop 138 END DO ! jj loop 139 140 CALL lbc_lnk( uwdlmt, 'U', 1. ) ; CALL lbc_lnk( vwdlmt , 'V', 1. ) 141 142 IF(zflag == 0) EXIT 143 115 zflag = 0 ! flag indicating if any further iteration is needed? 116 117 zflxu1(:,:) = zflxu(:,:) * wdlmt(:,:) 118 zflxv1(:,:) = zflxv(:,:) * wdlmt(:,:) 119 120 DO jj = 2, jpjm1 121 DO ji = fs_2, fs_jpim1 ! vector opt. 122 123 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 124 125 ztmp = e1t(ji,jj) * e2t(ji,jj) 126 127 zflxn(ji,jj) = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) + & 128 & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji, jj-1), 0._wp) 129 130 zdep1 = (zflxp(ji,jj) * wdlmt(ji,jj) + zflxn(ji,jj)) * z2dt / ztmp 131 zdep2 = bathy(ji,jj) + sshb(ji,jj) - rn_wadmin 132 133 IF(zdep1 > zdep2) THEN 134 zflag = 1 135 zcoef = ( ( zdep2 - rn_wadmine ) * ztmp - zflxn(ji,jj) * z2dt ) / ( zflxp(ji,jj) * z2dt ) 136 zcoef = max(zcoef, 0._wp) 137 IF(zflxu1(ji, jj ) >= 0._wp) wdlmt(ji, jj ) = zcoef 138 IF(zflxu1(ji-1,jj ) < 0._wp) wdlmt(ji-1,jj ) = zcoef 139 IF(zflxv1(ji, jj ) >= 0._wp) wdlmt(ji, jj ) = zcoef 140 IF(zflxv1(ji, jj-1) < 0._wp) wdlmt(ji, jj-1) = zcoef 141 END IF 142 END DO ! ji loop 143 END DO ! jj loop 144 145 CALL lbc_lnk( wdlmt, 'T', 1. ) 146 147 IF(zflag == 0) EXIT 148 149 IF(jk1 == nn_waditr) THEN 150 IF(zflxu1(ji, jj ) >= 0._wp) wdlmt(ji, jj ) = 0._wp 151 IF(zflxu1(ji-1,jj ) < 0._wp) wdlmt(ji-1,jj ) = 0._wp 152 IF(zflxv1(ji, jj ) >= 0._wp) wdlmt(ji, jj ) = 0._wp 153 IF(zflxv1(ji, jj-1) < 0._wp) wdlmt(ji, jj-1) = 0._wp 154 END IF 155 144 156 END DO ! jk1 loop 145 157 … … 147 159 DO jj = 2, jpjm1 148 160 DO ji = fs_2, fs_jpim1 ! vector opt. 149 un(ji, jj,jk) = 0.5_wp * ( sign(1._wp, un(ji,jj,jk)) + 1._wp ) * uwdlmt(ji,jj)150 vn(ji, jj,jk) = 0.5_wp * ( sign(1._wp, vn(ji,jj,jk)) + 1._wp ) * vwdlmt(ji,jj)151 un(ji-1,jj, jk) = 0.5_wp * ( 1._wp - sign(1._wp, un(ji-1,jj,jk) )) * uwdlmt(ji-1,jj)152 vn(ji, jj-1,jk) = 0.5_wp * ( 1._wp - sign(1._wp, vn(ji,jj-1,jk) )) * vwdlmt(ji,jj-1)161 un(ji, jj, jk) = ( sign(0.5_wp, un(ji, jj, jk)) + 0.5_wp ) * wdlmt(ji ,jj) 162 vn(ji, jj, jk) = ( sign(0.5_wp, vn(ji, jj, jk)) + 0.5_wp ) * wdlmt(ji ,jj) 163 un(ji-1,jj, jk) = (-sign(0.5_wp, un(ji-1,jj, jk)) + 0.5_wp ) * wdlmt(ji-1,jj) 164 vn(ji, jj-1,jk) = (-sign(0.5_wp, vn(ji, jj-1,jk)) + 0.5_wp ) * wdlmt(ji, jj-1) 153 165 END DO 154 166 END DO … … 161 173 ! 162 174 ! 163 CALL wrk_dealloc( jpi, jpj, zflx u, zflxv, zflxu1, zflxv1, uwdlmt, vwdlmt)175 CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1, wdlmt ) 164 176 ! 165 177 END IF
Note: See TracChangeset
for help on using the changeset viewer.
