- Timestamp:
- 2017-12-14T11:10:02+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90
r9019 r9031 17 17 18 18 PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90 19 19 #if defined key_vertical 20 PUBLIC reconstructandremap ! remapping routine 21 #endif 20 22 ! !!* Namelist namagrif: AGRIF parameters 21 23 LOGICAL , PUBLIC :: ln_spc_dyn = .FALSE. !: 22 INTEGER , PUBLIC :: nn_cln_update = 3 !: update frequency23 24 INTEGER , PUBLIC, PARAMETER :: nn_sponge_len = 2 !: Sponge width (in number of parent grid points) 24 25 REAL(wp), PUBLIC :: rn_sponge_tra = 2800. !: sponge coeff. for tracers 25 26 REAL(wp), PUBLIC :: rn_sponge_dyn = 2800. !: sponge coeff. for dynamics 26 27 LOGICAL , PUBLIC :: ln_chk_bathy = .FALSE. !: check of parent bathymetry 27 28 LOGICAL , PUBLIC :: lk_agrif_clp = .TRUE. !: Force clamped bcs 28 29 ! !!! OLD namelist names 29 INTEGER , PUBLIC :: nbcline = 0 !: update counter30 INTEGER , PUBLIC :: nbclineupdate !: update frequency31 30 REAL(wp), PUBLIC :: visc_tra !: sponge coeff. for tracers 32 31 REAL(wp), PUBLIC :: visc_dyn !: sponge coeff. for dynamics … … 35 34 LOGICAL , PUBLIC :: spongedoneU = .FALSE. !: dynamics sponge layer indicator 36 35 LOGICAL , PUBLIC :: lk_agrif_fstep = .TRUE. !: if true: first step 37 LOGICAL , PUBLIC :: lk_agrif_doupd = .TRUE. !: if true: send update from current grid38 36 LOGICAL , PUBLIC :: lk_agrif_debug = .FALSE. !: if true: print debugging info 39 37 … … 102 100 END FUNCTION agrif_oce_alloc 103 101 102 #if defined key_vertical 103 SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout) 104 !!---------------------------------------------------------------------- 105 !! *** FUNCTION reconstructandremap *** 106 !!---------------------------------------------------------------------- 107 IMPLICIT NONE 108 INTEGER N, Nout 109 REAL(wp) tabin(N), tabout(Nout) 110 REAL(wp) hin(N), hout(Nout) 111 REAL(wp) coeffremap(N,3),zwork(N,3) 112 REAL(wp) zwork2(N+1,3) 113 INTEGER jk 114 DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8 115 REAL(wp) q,q01,q02,q001,q002,q0 116 REAL(wp) z_win(1:N+1), z_wout(1:Nout+1) 117 REAL(wp),PARAMETER :: dpthin = 1.D-3 118 INTEGER :: k1, kbox, ktop, ka, kbot 119 REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 120 121 z_win(1)=0.; z_wout(1)= 0. 122 DO jk=1,N 123 z_win(jk+1)=z_win(jk)+hin(jk) 124 ENDDO 125 126 DO jk=1,Nout 127 z_wout(jk+1)=z_wout(jk)+hout(jk) 128 ENDDO 129 130 DO jk=2,N 131 zwork(jk,1)=1./(hin(jk-1)+hin(jk)) 132 ENDDO 133 134 DO jk=2,N-1 135 q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1)) 136 zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1) 137 zwork(jk,3)=q0 138 ENDDO 139 140 DO jk= 2,N 141 zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1)) 142 ENDDO 143 144 coeffremap(:,1) = tabin(:) 145 146 DO jk=2,N-1 147 q001 = hin(jk)*zwork2(jk+1,1) 148 q002 = hin(jk)*zwork2(jk,1) 149 IF (q001*q002 < 0) then 150 q001 = 0. 151 q002 = 0. 152 ENDIF 153 q=zwork(jk,2) 154 q01=q*zwork2(jk+1,1) 155 q02=q*zwork2(jk,1) 156 IF (abs(q001) > abs(q02)) q001 = q02 157 IF (abs(q002) > abs(q01)) q002 = q01 158 159 q=(q001-q002)*zwork(jk,3) 160 q001=q001-q*hin(jk+1) 161 q002=q002+q*hin(jk-1) 162 163 coeffremap(jk,3)=coeffremap(jk,1)+q001 164 coeffremap(jk,2)=coeffremap(jk,1)-q002 165 166 zwork2(jk,1)=(2.*q001-q002)**2 167 zwork2(jk,2)=(2.*q002-q001)**2 168 ENDDO 169 170 DO jk=1,N 171 IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN 172 coeffremap(jk,3) = coeffremap(jk,1) 173 coeffremap(jk,2) = coeffremap(jk,1) 174 zwork2(jk,1) = 0. 175 zwork2(jk,2) = 0. 176 ENDIF 177 ENDDO 178 179 DO jk=2,N 180 q002=max(zwork2(jk-1,2),dsmll) 181 q001=max(zwork2(jk,1),dsmll) 182 zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002) 183 ENDDO 184 185 zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3) 186 zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3) 187 188 DO jk=1,N 189 q01=zwork2(jk+1,3)-coeffremap(jk,1) 190 q02=coeffremap(jk,1)-zwork2(jk,3) 191 q001=2.*q01 192 q002=2.*q02 193 IF (q01*q02<0) then 194 q01=0. 195 q02=0. 196 ELSEIF (abs(q01)>abs(q002)) then 197 q01=q002 198 ELSEIF (abs(q02)>abs(q001)) then 199 q02=q001 200 ENDIF 201 coeffremap(jk,2)=coeffremap(jk,1)-q02 202 coeffremap(jk,3)=coeffremap(jk,1)+q01 203 ENDDO 204 205 zbot=0.0 206 kbot=1 207 DO jk=1,Nout 208 ztop=zbot !top is bottom of previous layer 209 ktop=kbot 210 IF (ztop.GE.z_win(ktop+1)) then 211 ktop=ktop+1 212 ENDIF 213 214 zbot=z_wout(jk+1) 215 zthk=zbot-ztop 216 217 IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN 218 219 kbot=ktop 220 DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) 221 kbot=kbot+1 222 ENDDO 223 zbox=zbot 224 DO k1= jk+1,Nout 225 IF (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN 226 exit !thick layer 227 ELSE 228 zbox=z_wout(k1+1) !include thin adjacent layers 229 IF(zbox.EQ.z_wout(Nout+1)) THEN 230 exit !at bottom 231 ENDIF 232 ENDIF 233 ENDDO 234 zthk=zbox-ztop 235 236 kbox=ktop 237 DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) 238 kbox=kbox+1 239 ENDDO 240 241 IF(ktop.EQ.kbox) THEN 242 IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN 243 IF(hin(kbox).GT.dpthin) THEN 244 q001 = (zbox-z_win(kbox))/hin(kbox) 245 q002 = (ztop-z_win(kbox))/hin(kbox) 246 q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) 247 q02=q01-1.+(q001+q002) 248 q0=1.-q01-q02 249 ELSE 250 q0 = 1.0 251 q01 = 0. 252 q02 = 0. 253 ENDIF 254 tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 255 ELSE 256 tabout(jk) = tabin(kbox) 257 ENDIF 258 ELSE 259 IF(ktop.LE.jk .AND. kbox.GE.jk) THEN 260 ka = jk 261 ELSEIF (kbox-ktop.GE.3) THEN 262 ka = (kbox+ktop)/2 263 ELSEIF (hin(ktop).GE.hin(kbox)) THEN 264 ka = ktop 265 ELSE 266 ka = kbox 267 ENDIF !choose ka 268 269 offset=coeffremap(ka,1) 270 271 qtop = z_win(ktop+1)-ztop !partial layer thickness 272 IF(hin(ktop).GT.dpthin) THEN 273 q=(ztop-z_win(ktop))/hin(ktop) 274 q01=q*(q-1.) 275 q02=q01+q 276 q0=1-q01-q02 277 ELSE 278 q0 = 1. 279 q01 = 0. 280 q02 = 0. 281 ENDIF 282 283 tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop 284 285 DO k1= ktop+1,kbox-1 286 tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) 287 ENDDO !k1 288 289 qbot = zbox-z_win(kbox) !partial layer thickness 290 IF(hin(kbox).GT.dpthin) THEN 291 q=qbot/hin(kbox) 292 q01=(q-1.)**2 293 q02=q01-1.+q 294 q0=1-q01-q02 295 ELSE 296 q0 = 1.0 297 q01 = 0. 298 q02 = 0. 299 ENDIF 300 301 tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot 302 303 rpsum=1.0d0/zthk 304 tabout(jk)=offset+tsum*rpsum 305 306 ENDIF !single or multiple layers 307 ELSE 308 IF (jk==1) THEN 309 write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1) 310 ENDIF 311 tabout(jk) = tabout(jk-1) 312 313 ENDIF !normal:thin layer 314 ENDDO !jk 315 316 return 317 end subroutine reconstructandremap 318 #endif 319 104 320 #endif 105 321 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
