Changeset 11590 for NEMO/branches/2019/dev_r11233_AGRIF05_jchanut_vert_coord_interp/src/NST/agrif_oce.F90
 Timestamp:
 20190923T18:25:29+02:00 (14 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

NEMO/branches/2019/dev_r11233_AGRIF05_jchanut_vert_coord_interp/src/NST/agrif_oce.F90
r11574 r11590 17 17 18 18 PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90 19 #if defined key_vertical 20 PUBLIC reconstructandremap ! remapping routine 21 #endif 19 22 20 ! !!* Namelist namagrif: AGRIF parameters 23 21 LOGICAL , PUBLIC :: ln_agrif_2way = .TRUE. !: activate two way nesting … … 95 93 END FUNCTION agrif_oce_alloc 96 94 97 #if defined key_vertical98 SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)99 !!100 !! *** FUNCTION reconstructandremap ***101 !!102 IMPLICIT NONE103 INTEGER N, Nout104 REAL(wp) tabin(N), tabout(Nout)105 REAL(wp) hin(N), hout(Nout)106 REAL(wp) coeffremap(N,3),zwork(N,3)107 REAL(wp) zwork2(N+1,3)108 INTEGER jk109 DOUBLE PRECISION, PARAMETER :: dsmll=1.0d8110 REAL(wp) q,q01,q02,q001,q002,q0111 REAL(wp) z_win(1:N+1), z_wout(1:Nout+1)112 REAL(wp),PARAMETER :: dpthin = 1.D3113 INTEGER :: k1, kbox, ktop, ka, kbot114 REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop115 116 z_win(1)=0.; z_wout(1)= 0.117 DO jk=1,N118 z_win(jk+1)=z_win(jk)+hin(jk)119 ENDDO120 121 DO jk=1,Nout122 z_wout(jk+1)=z_wout(jk)+hout(jk)123 ENDDO124 125 DO jk=2,N126 zwork(jk,1)=1./(hin(jk1)+hin(jk))127 ENDDO128 129 DO jk=2,N1130 q0 = 1./(hin(jk1)+hin(jk)+hin(jk+1))131 zwork(jk,2)=hin(jk1)+2.*hin(jk)+hin(jk+1)132 zwork(jk,3)=q0133 ENDDO134 135 DO jk= 2,N136 zwork2(jk,1)=zwork(jk,1)*(tabin(jk)tabin(jk1))137 ENDDO138 139 coeffremap(:,1) = tabin(:)140 141 DO jk=2,N1142 q001 = hin(jk)*zwork2(jk+1,1)143 q002 = hin(jk)*zwork2(jk,1)144 IF (q001*q002 < 0) then145 q001 = 0.146 q002 = 0.147 ENDIF148 q=zwork(jk,2)149 q01=q*zwork2(jk+1,1)150 q02=q*zwork2(jk,1)151 IF (abs(q001) > abs(q02)) q001 = q02152 IF (abs(q002) > abs(q01)) q002 = q01153 154 q=(q001q002)*zwork(jk,3)155 q001=q001q*hin(jk+1)156 q002=q002+q*hin(jk1)157 158 coeffremap(jk,3)=coeffremap(jk,1)+q001159 coeffremap(jk,2)=coeffremap(jk,1)q002160 161 zwork2(jk,1)=(2.*q001q002)**2162 zwork2(jk,2)=(2.*q002q001)**2163 ENDDO164 165 DO jk=1,N166 IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN167 coeffremap(jk,3) = coeffremap(jk,1)168 coeffremap(jk,2) = coeffremap(jk,1)169 zwork2(jk,1) = 0.170 zwork2(jk,2) = 0.171 ENDIF172 ENDDO173 174 DO jk=2,N175 q002=max(zwork2(jk1,2),dsmll)176 q001=max(zwork2(jk,1),dsmll)177 zwork2(jk,3)=(q001*coeffremap(jk1,3)+q002*coeffremap(jk,2))/(q001+q002)178 ENDDO179 180 zwork2(1,3) = 2*coeffremap(1,1)zwork2(2,3)181 zwork2(N+1,3)=2*coeffremap(N,1)zwork2(N,3)182 183 DO jk=1,N184 q01=zwork2(jk+1,3)coeffremap(jk,1)185 q02=coeffremap(jk,1)zwork2(jk,3)186 q001=2.*q01187 q002=2.*q02188 IF (q01*q02<0) then189 q01=0.190 q02=0.191 ELSEIF (abs(q01)>abs(q002)) then192 q01=q002193 ELSEIF (abs(q02)>abs(q001)) then194 q02=q001195 ENDIF196 coeffremap(jk,2)=coeffremap(jk,1)q02197 coeffremap(jk,3)=coeffremap(jk,1)+q01198 ENDDO199 200 zbot=0.0201 kbot=1202 DO jk=1,Nout203 ztop=zbot !top is bottom of previous layer204 ktop=kbot205 IF (ztop.GE.z_win(ktop+1)) then206 ktop=ktop+1207 ENDIF208 209 zbot=z_wout(jk+1)210 zthk=zbotztop211 212 IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN213 214 kbot=ktop215 DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)216 kbot=kbot+1217 ENDDO218 zbox=zbot219 DO k1= jk+1,Nout220 IF (z_wout(k1+1)z_wout(k1).GT.dpthin) THEN221 exit !thick layer222 ELSE223 zbox=z_wout(k1+1) !include thin adjacent layers224 IF(zbox.EQ.z_wout(Nout+1)) THEN225 exit !at bottom226 ENDIF227 ENDIF228 ENDDO229 zthk=zboxztop230 231 kbox=ktop232 DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)233 kbox=kbox+1234 ENDDO235 236 IF(ktop.EQ.kbox) THEN237 IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN238 IF(hin(kbox).GT.dpthin) THEN239 q001 = (zboxz_win(kbox))/hin(kbox)240 q002 = (ztopz_win(kbox))/hin(kbox)241 q01=q001**2+q002**2+q001*q002+1.2.*(q001+q002)242 q02=q011.+(q001+q002)243 q0=1.q01q02244 ELSE245 q0 = 1.0246 q01 = 0.247 q02 = 0.248 ENDIF249 tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)250 ELSE251 tabout(jk) = tabin(kbox)252 ENDIF253 ELSE254 IF(ktop.LE.jk .AND. kbox.GE.jk) THEN255 ka = jk256 ELSEIF (kboxktop.GE.3) THEN257 ka = (kbox+ktop)/2258 ELSEIF (hin(ktop).GE.hin(kbox)) THEN259 ka = ktop260 ELSE261 ka = kbox262 ENDIF !choose ka263 264 offset=coeffremap(ka,1)265 266 qtop = z_win(ktop+1)ztop !partial layer thickness267 IF(hin(ktop).GT.dpthin) THEN268 q=(ztopz_win(ktop))/hin(ktop)269 q01=q*(q1.)270 q02=q01+q271 q0=1q01q02272 ELSE273 q0 = 1.274 q01 = 0.275 q02 = 0.276 ENDIF277 278 tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))offset)*qtop279 280 DO k1= ktop+1,kbox1281 tsum =tsum +(coeffremap(k1,1)offset)*hin(k1)282 ENDDO !k1283 284 qbot = zboxz_win(kbox) !partial layer thickness285 IF(hin(kbox).GT.dpthin) THEN286 q=qbot/hin(kbox)287 q01=(q1.)**2288 q02=q011.+q289 q0=1q01q02290 ELSE291 q0 = 1.0292 q01 = 0.293 q02 = 0.294 ENDIF295 296 tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))offset)*qbot297 298 rpsum=1.0d0/zthk299 tabout(jk)=offset+tsum*rpsum300 301 ENDIF !single or multiple layers302 ELSE303 IF (jk==1) THEN304 write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1)305 ENDIF306 tabout(jk) = tabout(jk1)307 308 ENDIF !normal:thin layer309 ENDDO !jk310 311 return312 end subroutine reconstructandremap313 #endif314 315 95 #endif 316 96 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.