[Dart-dev] [3309] DART/trunk/models/wrf/model_mod.f90: Adding Ryan
Torn' s code for a forward operator for
nancy at subversion.ucar.edu
nancy at subversion.ucar.edu
Fri Apr 11 14:53:43 MDT 2008
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20080411/0d7460ec/attachment.html
-------------- next part --------------
Modified: DART/trunk/models/wrf/model_mod.f90
===================================================================
--- DART/trunk/models/wrf/model_mod.f90 2008-04-10 20:30:33 UTC (rev 3308)
+++ DART/trunk/models/wrf/model_mod.f90 2008-04-11 20:53:42 UTC (rev 3309)
@@ -1175,7 +1175,7 @@
real(r8), allocatable, dimension(:,:) :: psea, pp, pd
real(r8), allocatable, dimension(:) :: x1d,y1d,xx1d,yy1d
integer :: xlen, ylen, xxlen, yylen, ii1, ii2
-real(r8) :: clat, clon, cxmin, cymin
+real(r8) :: clat, clon, cxmin, cymin, maxwspd
! center_track_*** used to define center search area
integer :: center_track_xmin, center_track_ymin, &
@@ -2255,74 +2255,95 @@
allocate(y1d(ylen))
allocate(xx1d(xxlen))
allocate(yy1d(yylen))
+ maxwspd = missing_r8
-!! compute sea-level pressure
- do i1 = center_track_xmin, center_track_xmax
- do i2 = center_track_ymin, center_track_ymax
- do k2 = 1,wrf%dom(id)%var_size(3,TYPE_T)
- p1d(k2) = model_pressure_t(i1,i2,k2,id,x)
- t1d(k2) = ts0 + x(wrf%dom(id)%dart_ind(i1,i2,k2,TYPE_T))
- qv1d(k2)= x(wrf%dom(id)%dart_ind(i1,i2,k2,TYPE_QV))
- z1d(k2) = ( x(wrf%dom(id)%dart_ind(i1,i2,k2,TYPE_GZ))+wrf%dom(id)%phb(i1,i2,k2) + &
- x(wrf%dom(id)%dart_ind(i1,i2,k2+1,TYPE_GZ))+wrf%dom(id)%phb(i1,i2,k2+1) &
- )*0.5_r8/gravity
- enddo
- call compute_seaprs(wrf%dom(id)%bt, z1d, t1d, p1d, qv1d, &
- psea(i1-center_track_xmin+1,i2-center_track_ymin+1),debug)
- enddo
- enddo
+!! compute sea-level pressure and max wind speed
+ if ( obs_kind == KIND_VORTEX_WMAX ) then
-!! spline-interpolation
- do i1 = 1,xlen
- x1d(i1) = (i1-1)+center_track_xmin
- enddo
- do i2 = 1,ylen
- y1d(i2) = (i2-1)+center_track_ymin
- enddo
- do ii1 = 1,xxlen
- xx1d(ii1) = center_track_xmin+real(ii1-1,r8)*1_r8/real(center_spline_grid_scale,r8)
- enddo
- do ii2 = 1,yylen
- yy1d(ii2) = center_track_ymin+real(ii2-1,r8)*1_r8/real(center_spline_grid_scale,r8)
- enddo
+ do i1 = center_track_xmin, center_track_xmax
+ do i2 = center_track_ymin, center_track_ymax
- call splie2(x1d,y1d,psea,xlen,ylen,pd)
+ if ( wrf%dom(id)%surf_obs ) then
+ maxwspd = max(maxwspd, sqrt( x(wrf%dom(id)%dart_ind(i1,i2,1,TYPE_U10)) ** 2 + &
+ x(wrf%dom(id)%dart_ind(i1,i2,1,TYPE_V10)) ** 2 ))
+ else
+ maxwspd = max(maxwspd,sqrt((0.5_r8*(x(wrf%dom(id)%dart_ind(i1, i2,1,TYPE_U)) + &
+ x(wrf%dom(id)%dart_ind(i1+1,i2,1,TYPE_U))))**2 + &
+ (0.5_r8*(x(wrf%dom(id)%dart_ind(i1,i2, 1,TYPE_V)) + &
+ x(wrf%dom(id)%dart_ind(i1,i2+1,1,TYPE_V))))**2))
+ end if
- pres1 = 1.e20
- cxmin = -1
- cymin = -1
- do ii1=1,xxlen
- do ii2=1,yylen
- call splin2(x1d,y1d,psea,pd,xlen,ylen,xx1d(ii1),yy1d(ii2),pp(ii1,ii2))
- if (pres1 .gt. pp(ii1,ii2)) then
- pres1=pp(ii1,ii2)
- cxmin=xx1d(ii1)
- cymin=yy1d(ii2)
- endif
- enddo
- enddo
+ enddo
+ enddo
+ fld(1) = maxwspd
-!! if too close to the edge of the box, reset to observed center
- if( cxmin-xx1d(1) < 1_r8 .or. xx1d(xxlen)-cxmin < 1_r8 .or. &
- cymin-yy1d(1) < 1_r8 .or. yy1d(yylen)-cymin < 1_r8 ) then
- cxmin = xloc
- cymin = yloc
- call splin2(x1d,y1d,psea,pd,xlen,ylen,cxmin,cymin,pres1)
- endif
+ else
- call ij_to_latlon(wrf%dom(id)%proj, cxmin, cymin, clat, clon)
+ do i1 = center_track_xmin, center_track_xmax
+ do i2 = center_track_ymin, center_track_ymax
+ do k2 = 1,wrf%dom(id)%var_size(3,TYPE_T)
+ p1d(k2) = model_pressure_t(i1,i2,k2,id,x)
+ t1d(k2) = ts0 + x(wrf%dom(id)%dart_ind(i1,i2,k2,TYPE_T))
+ qv1d(k2)= x(wrf%dom(id)%dart_ind(i1,i2,k2,TYPE_QV))
+ z1d(k2) = ( x(wrf%dom(id)%dart_ind(i1,i2,k2,TYPE_GZ))+wrf%dom(id)%phb(i1,i2,k2) + &
+ x(wrf%dom(id)%dart_ind(i1,i2,k2+1,TYPE_GZ))+wrf%dom(id)%phb(i1,i2,k2+1) &
+ )*0.5_r8/gravity
+ enddo
+ call compute_seaprs(wrf%dom(id)%bt, z1d, t1d, p1d, qv1d, &
+ psea(i1-center_track_xmin+1,i2-center_track_ymin+1),debug)
+ enddo
+ enddo
- if( obs_kind == KIND_VORTEX_LAT ) then
- fld(1) = clat
- else if( obs_kind == KIND_VORTEX_LON ) then
- fld(1) = clon
- else if( obs_kind == KIND_VORTEX_PMIN ) then
- fld(1) = pres1
- else if( obs_kind == KIND_VORTEX_WMAX ) then
- fld(1) = missing_r8
- ! not implemented yet
- endif
+!! spline-interpolation
+ do i1 = 1,xlen
+ x1d(i1) = (i1-1)+center_track_xmin
+ enddo
+ do i2 = 1,ylen
+ y1d(i2) = (i2-1)+center_track_ymin
+ enddo
+ do ii1 = 1,xxlen
+ xx1d(ii1) = center_track_xmin+real(ii1-1,r8)*1_r8/real(center_spline_grid_scale,r8)
+ enddo
+ do ii2 = 1,yylen
+ yy1d(ii2) = center_track_ymin+real(ii2-1,r8)*1_r8/real(center_spline_grid_scale,r8)
+ enddo
+ call splie2(x1d,y1d,psea,xlen,ylen,pd)
+
+ pres1 = 1.e20
+ cxmin = -1
+ cymin = -1
+ do ii1=1,xxlen
+ do ii2=1,yylen
+ call splin2(x1d,y1d,psea,pd,xlen,ylen,xx1d(ii1),yy1d(ii2),pp(ii1,ii2))
+ if (pres1 .gt. pp(ii1,ii2)) then
+ pres1=pp(ii1,ii2)
+ cxmin=xx1d(ii1)
+ cymin=yy1d(ii2)
+ endif
+ enddo
+ enddo
+
+!! if too close to the edge of the box, reset to observed center
+ if( cxmin-xx1d(1) < 1_r8 .or. xx1d(xxlen)-cxmin < 1_r8 .or. &
+ cymin-yy1d(1) < 1_r8 .or. yy1d(yylen)-cymin < 1_r8 ) then
+ cxmin = xloc
+ cymin = yloc
+ call splin2(x1d,y1d,psea,pd,xlen,ylen,cxmin,cymin,pres1)
+ endif
+
+ call ij_to_latlon(wrf%dom(id)%proj, cxmin, cymin, clat, clon)
+
+ if( obs_kind == KIND_VORTEX_LAT ) then
+ fld(1) = clat
+ else if( obs_kind == KIND_VORTEX_LON ) then
+ fld(1) = clon
+ else if( obs_kind == KIND_VORTEX_PMIN ) then
+ fld(1) = pres1
+ endif
+
+ endif
+
deallocate(p1d, t1d, qv1d, z1d)
deallocate(psea,pd,pp,x1d,y1d,xx1d,yy1d)
More information about the Dart-dev
mailing list