Changeset 8855
 Timestamp:
 20171130T12:53:57+01:00 (7 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90
r8135 r8855 105 105 END FUNCTION agrif_oce_alloc 106 106 107 subroutine reconstructandremap(tabin,hin,tabout,hout,N,Nout) 108 implicit none 109 integer N, Nout 110 real tabin(N), tabout(Nout) 111 real hin(N), hout(Nout) 112 real coeffremap(N,3),zwork(N,3) 113 real zwork2(N+1,3) 114 integer k 115 double precision, parameter :: dsmll=1.0d8 116 real q,q01,q02,q001,q002,q0 117 real z_win(1:N+1), z_wout(1:Nout+1) 118 real,parameter :: dpthin = 1.D3 119 integer :: k1, kbox, ktop, ka, kbot 120 real :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 107 SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout) 108 !! 109 !! *** FUNCTION reconstructandremap *** 110 !! 111 IMPLICIT NONE 112 INTEGER N, Nout 113 REAL(wp) tabin(N), tabout(Nout) 114 REAL(wp) hin(N), hout(Nout) 115 REAL(wp) coeffremap(N,3),zwork(N,3) 116 REAL(wp) zwork2(N+1,3) 117 INTEGER jk 118 DOUBLE PRECISION, PARAMETER :: dsmll=1.0d8 119 REAL(wp) q,q01,q02,q001,q002,q0 120 REAL(wp) z_win(1:N+1), z_wout(1:Nout+1) 121 REAL(wp),PARAMETER :: dpthin = 1.D3 122 INTEGER :: k1, kbox, ktop, ka, kbot 123 REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 121 124 122 125 z_win(1)=0.; z_wout(1)= 0. 123 dok=1,N124 z_win(k+1)=z_win(k)+hin(k)125 enddo126 DO jk=1,N 127 z_win(jk+1)=z_win(jk)+hin(jk) 128 ENDDO 126 129 127 dok=1,Nout128 z_wout(k+1)=z_wout(k)+hout(k)129 enddo130 131 dok=2,N132 zwork(k,1)=1./(hin(k1)+hin(k))133 enddo134 135 dok=2,N1136 q0 = 1./(hin(k1)+hin(k)+hin(k+1))137 zwork(k,2)=hin(k1)+2.*hin(k)+hin(k+1)138 zwork(k,3)=q0139 enddo140 141 dok= 2,N142 zwork2(k,1)=zwork(k,1)*(tabin(k)tabin(k1))143 enddo144 145 130 DO jk=1,Nout 131 z_wout(jk+1)=z_wout(jk)+hout(jk) 132 ENDDO 133 134 DO jk=2,N 135 zwork(jk,1)=1./(hin(jk1)+hin(jk)) 136 ENDDO 137 138 DO jk=2,N1 139 q0 = 1./(hin(jk1)+hin(jk)+hin(jk+1)) 140 zwork(jk,2)=hin(jk1)+2.*hin(jk)+hin(jk+1) 141 zwork(jk,3)=q0 142 ENDDO 143 144 DO jk= 2,N 145 zwork2(jk,1)=zwork(jk,1)*(tabin(jk)tabin(jk1)) 146 ENDDO 147 148 coeffremap(:,1) = tabin(:) 146 149 147 dok=2,N1148 q001 = hin(k)*zwork2(k+1,1)149 q002 = hin(k)*zwork2(k,1)150 if(q001*q002 < 0) then151 q001 = 0.152 q002 = 0.153 endif154 q=zwork(k,2)155 q01=q*zwork2(k+1,1)156 q02=q*zwork2(k,1)157 if(abs(q001) > abs(q02)) q001 = q02158 if(abs(q002) > abs(q01)) q002 = q01159 160 q=(q001q002)*zwork(k,3)161 q001=q001q*hin(k+1)162 q002=q002+q*hin(k1)163 164 coeffremap(k,3)=coeffremap(k,1)+q001165 coeffremap(k,2)=coeffremap(k,1)q002166 167 zwork2(k,1)=(2.*q001q002)**2168 zwork2(k,2)=(2.*q002q001)**2169 enddo170 171 dok=1,N172 if (k.eq.1 .or. k.eq.N .or. hin(k).le.dpthin) then173 coeffremap(k,3) = coeffremap(k,1)174 coeffremap(k,2) = coeffremap(k,1)175 zwork2(k,1) = 0.176 zwork2(k,2) = 0.177 endif178 enddo179 180 dok=2,N181 q002=max(zwork2(k1,2),dsmll)182 q001=max(zwork2(k,1),dsmll)183 zwork2(k,3)=(q001*coeffremap(k1,3)+q002*coeffremap(k,2))/(q001+q002)184 enddo185 186 187 150 DO jk=2,N1 151 q001 = hin(jk)*zwork2(jk+1,1) 152 q002 = hin(jk)*zwork2(jk,1) 153 IF (q001*q002 < 0) then 154 q001 = 0. 155 q002 = 0. 156 ENDIF 157 q=zwork(jk,2) 158 q01=q*zwork2(jk+1,1) 159 q02=q*zwork2(jk,1) 160 IF (abs(q001) > abs(q02)) q001 = q02 161 IF (abs(q002) > abs(q01)) q002 = q01 162 163 q=(q001q002)*zwork(jk,3) 164 q001=q001q*hin(jk+1) 165 q002=q002+q*hin(jk1) 166 167 coeffremap(jk,3)=coeffremap(jk,1)+q001 168 coeffremap(jk,2)=coeffremap(jk,1)q002 169 170 zwork2(jk,1)=(2.*q001q002)**2 171 zwork2(jk,2)=(2.*q002q001)**2 172 ENDDO 173 174 DO jk=1,N 175 IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN 176 coeffremap(jk,3) = coeffremap(jk,1) 177 coeffremap(jk,2) = coeffremap(jk,1) 178 zwork2(jk,1) = 0. 179 zwork2(jk,2) = 0. 180 ENDIF 181 ENDDO 182 183 DO jk=2,N 184 q002=max(zwork2(jk1,2),dsmll) 185 q001=max(zwork2(jk,1),dsmll) 186 zwork2(jk,3)=(q001*coeffremap(jk1,3)+q002*coeffremap(jk,2))/(q001+q002) 187 ENDDO 188 189 zwork2(1,3) = 2*coeffremap(1,1)zwork2(2,3) 190 zwork2(N+1,3)=2*coeffremap(N,1)zwork2(N,3) 188 191 189 dok=1,N190 q01=zwork2(k+1,3)coeffremap(k,1)191 q02=coeffremap(k,1)zwork2(k,3)192 q001=2.*q01193 q002=2.*q02194 if(q01*q02<0) then195 q01=0.196 q02=0.197 elseif(abs(q01)>abs(q002)) then198 q01=q002199 elseif(abs(q02)>abs(q001)) then200 q02=q001201 endif202 coeffremap(k,2)=coeffremap(k,1)q02203 coeffremap(k,3)=coeffremap(k,1)+q01204 enddo192 DO jk=1,N 193 q01=zwork2(jk+1,3)coeffremap(jk,1) 194 q02=coeffremap(jk,1)zwork2(jk,3) 195 q001=2.*q01 196 q002=2.*q02 197 IF (q01*q02<0) then 198 q01=0. 199 q02=0. 200 ELSEIF (abs(q01)>abs(q002)) then 201 q01=q002 202 ELSEIF (abs(q02)>abs(q001)) then 203 q02=q001 204 ENDIF 205 coeffremap(jk,2)=coeffremap(jk,1)q02 206 coeffremap(jk,3)=coeffremap(jk,1)+q01 207 ENDDO 205 208 206 209 zbot=0.0 207 210 kbot=1 208 dok=1,Nout209 ztop=zbot !top is bottom of previous layer210 ktop=kbot211 if (ztop.ge.z_win(ktop+1)) then212 ktop=ktop+1213 endif214 215 zbot=z_wout(k+1)216 zthk=zbotztop217 218 if (zthk.gt.dpthin .and. ztop.lt.z_wout(Nout+1)) then219 220 kbot=ktop221 dowhile (z_win(kbot+1).lt.zbot.and.kbot.lt.N)222 kbot=kbot+1223 enddo224 zbox=zbot225 do k1=k+1,Nout226 if (z_wout(k1+1)z_wout(k1).gt.dpthin) then227 exit !thick layer228 else229 zbox=z_wout(k1+1) !include thin adjacent layers230 if (zbox.eq.z_wout(Nout+1)) then231 exit !at bottom232 endif233 endif234 enddo235 zthk=zboxztop236 237 kbox=ktop238 dowhile (z_win(kbox+1).lt.zbox.and.kbox.lt.N)239 kbox=kbox+1240 enddo211 DO jk=1,Nout 212 ztop=zbot !top is bottom of previous layer 213 ktop=kbot 214 IF (ztop.GE.z_win(ktop+1)) then 215 ktop=ktop+1 216 ENDIF 217 218 zbot=z_wout(jk+1) 219 zthk=zbotztop 220 221 IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN 222 223 kbot=ktop 224 DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) 225 kbot=kbot+1 226 ENDDO 227 zbox=zbot 228 DO k1= jk+1,Nout 229 IF (z_wout(k1+1)z_wout(k1).GT.dpthin) THEN 230 exit !thick layer 231 ELSE 232 zbox=z_wout(k1+1) !include thin adjacent layers 233 IF(zbox.EQ.z_wout(Nout+1)) THEN 234 exit !at bottom 235 ENDIF 236 ENDIF 237 ENDDO 238 zthk=zboxztop 239 240 kbox=ktop 241 DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) 242 kbox=kbox+1 243 ENDDO 241 244 242 if (ktop.eq.kbox) then 243 244 245 if (z_wout(k) .ne.z_win(kbox) .or.z_wout(k+1).ne.z_win(kbox+1) ) then 246 247 if (hin(kbox).gt.dpthin) then 248 q001 = (zboxz_win(kbox))/hin(kbox) 249 q002 = (ztopz_win(kbox))/hin(kbox) 250 q01=q001**2+q002**2+q001*q002+1.2.*(q001+q002) 251 q02=q011.+(q001+q002) 252 q0=1.q01q02 253 else 254 q0 = 1.0 255 q01 = 0. 256 q02 = 0. 257 endif 258 tabout(k)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 245 IF(ktop.EQ.kbox) THEN 246 IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN 247 IF(hin(kbox).GT.dpthin) THEN 248 q001 = (zboxz_win(kbox))/hin(kbox) 249 q002 = (ztopz_win(kbox))/hin(kbox) 250 q01=q001**2+q002**2+q001*q002+1.2.*(q001+q002) 251 q02=q011.+(q001+q002) 252 q0=1.q01q02 253 ELSE 254 q0 = 1.0 255 q01 = 0. 256 q02 = 0. 257 ENDIF 258 tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 259 ELSE 260 tabout(jk) = tabin(kbox) 261 ENDIF 262 ELSE 263 IF(ktop.LE.jk .AND. kbox.GE.jk) THEN 264 ka = jk 265 ELSEIF (kboxktop.GE.3) THEN 266 ka = (kbox+ktop)/2 267 ELSEIF (hin(ktop).GE.hin(kbox)) THEN 268 ka = ktop 269 ELSE 270 ka = kbox 271 ENDIF !choose ka 272 273 offset=coeffremap(ka,1) 274 275 qtop = z_win(ktop+1)ztop !partial layer thickness 276 IF(hin(ktop).GT.dpthin) THEN 277 q=(ztopz_win(ktop))/hin(ktop) 278 q01=q*(q1.) 279 q02=q01+q 280 q0=1q01q02 281 ELSE 282 q0 = 1. 283 q01 = 0. 284 q02 = 0. 285 ENDIF 286 287 tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))offset)*qtop 288 289 DO k1= ktop+1,kbox1 290 tsum =tsum +(coeffremap(k1,1)offset)*hin(k1) 291 ENDDO !k1 292 293 qbot = zboxz_win(kbox) !partial layer thickness 294 IF(hin(kbox).GT.dpthin) THEN 295 q=qbot/hin(kbox) 296 q01=(q1.)**2 297 q02=q011.+q 298 q0=1q01q02 299 ELSE 300 q0 = 1.0 301 q01 = 0. 302 q02 = 0. 303 ENDIF 259 304 260 else 261 tabout(k) = tabin(kbox) 262 263 endif 264 265 else 266 267 if (ktop.le.k .and. kbox.ge.k) then 268 ka = k 269 elseif (kboxktop.ge.3) then 270 ka = (kbox+ktop)/2 271 elseif (hin(ktop).ge.hin(kbox)) then 272 ka = ktop 273 else 274 ka = kbox 275 endif !choose ka 276 277 offset=coeffremap(ka,1) 278 279 qtop = z_win(ktop+1)ztop !partial layer thickness 280 if (hin(ktop).gt.dpthin) then 281 q=(ztopz_win(ktop))/hin(ktop) 282 q01=q*(q1.) 283 q02=q01+q 284 q0=1q01q02 285 else 286 q0 = 1. 287 q01 = 0. 288 q02 = 0. 289 endif 290 291 tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))offset)*qtop 292 293 do k1= ktop+1,kbox1 294 tsum =tsum +(coeffremap(k1,1)offset)*hin(k1) 295 enddo !k1 296 297 298 qbot = zboxz_win(kbox) !partial layer thickness 299 if (hin(kbox).gt.dpthin) then 300 q=qbot/hin(kbox) 301 q01=(q1.)**2 302 q02=q011.+q 303 q0=1q01q02 304 else 305 q0 = 1.0 306 q01 = 0. 307 q02 = 0. 308 endif 309 310 tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))offset)*qbot 311 312 rpsum=1.0d0/zthk 313 tabout(k)=offset+tsum*rpsum 314 315 endif !single or multiple layers 316 else 317 if (k==1) then 318 write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(k+1),hout(1) 319 endif 320 tabout(k) = tabout(k1) 321 322 endif !normal:thin layer 323 enddo !k 305 tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))offset)*qbot 306 307 rpsum=1.0d0/zthk 308 tabout(jk)=offset+tsum*rpsum 309 310 ENDIF !single or multiple layers 311 ELSE 312 IF (jk==1) THEN 313 write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1) 314 ENDIF 315 tabout(jk) = tabout(jk1) 316 317 ENDIF !normal:thin layer 318 ENDDO !jk 324 319 325 320 return
Note: See TracChangeset
for help on using the changeset viewer.