MODULE v5d_module USE gridinfo_module USE output_module USE map_utils USE misc_definitions_module CONTAINS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: v5d_varnames ! Purpose: Get info about the variables to send to v5d file ! This routine assumes that the user only asked for variables ! that are available. If variables are not available v5d will no ! be working correctly. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE v5d_varnames (v5d_nl) IMPLICIT NONE ! ------ some vis5d stuff, contain MAXVARS,MAXTIMES,MAXROWS,MAXCOLUMNS,MAXLEVELS include 'v5df.h' integer, dimension(MAXLEVELS) :: v5d_nl character (len=4000) :: v5d_fields integer :: is_there, is_3d character (len=1000) :: known_3d_fields character (len=20) :: dummy known_3d_fields = 'P,PB,QCLOUD,QGRAUP,QICE,QNICE,QRAIN,QSNOW, & &QVAPOR,T,TKE,TKE_MYJ,U,V,PH,PHB,W,cape,cin,dbz,height,pressure,rh, & &tc,td,theta,tk,umet,vmet,wdir,wspd' varnames = ' ' numvars = 0 v5d_fields = plot_these_fields(2:len_trim(plot_these_fields)) is_there = INDEX(v5d_fields,",") DO WHILE ( is_there > 1 ) numvars = numvars + 1 varnames(numvars) = v5d_fields(1:is_there-1) dummy = ","//trim(varnames(numvars))//"," IF ( trim(varnames(numvars)) == 'SH2O' .OR. & trim(varnames(numvars)) == 'SMOIS' .OR. & trim(varnames(numvars)) == 'TSLB' ) THEN print*," WARING: We cannot deal with field ", trim(varnames(numvars)) print*," It will be removed from the plot list" is_3d = INDEX(plot_these_fields,trim(dummy)) plot_these_fields = plot_these_fields(1:is_3d)//plot_these_fields(is_3d+len_trim(dummy):len_trim(plot_these_fields)) varnames(numvars) = ' ' numvars = numvars - 1 ELSE v5d_nl(numvars) = 1 is_3d = INDEX(known_3d_fields,trim(varnames(numvars))) IF ( is_3d > 0 ) v5d_nl(numvars) = number_of_zlevs END IF v5d_fields = v5d_fields(is_there+1:len_trim(v5d_fields)) !!! Check for doubles is_there = INDEX(v5d_fields,trim(dummy)) DO WHILE ( is_there /= 0 ) v5d_fields = v5d_fields(1:is_there)//v5d_fields(is_there+len_trim(dummy):len_trim(v5d_fields)) is_there = INDEX(v5d_fields,trim(dummy)) END DO is_there = INDEX(v5d_fields,",") END DO END SUBROUTINE v5d_varnames !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: v5d_times ! Purpose: Get info about the variables to send to v5d file ! This routine assumes that the user only asked for variables ! that are available. If variables are not available v5d will no ! be working correctly. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE v5d_times (n_times, timestamp, datestamp) USE date_pack IMPLICIT NONE ! ------ some vis5d stuff, contain MAXVARS,MAXTIMES,MAXROWS,MAXCOLUMNS,MAXLEVELS include 'v5df.h' integer, dimension(MAXTIMES) :: timestamp, datestamp character (len=19) :: next_date integer :: tt, n_times integer :: year, month, day, hours, minutes, seconds datestamp = IMISSING timestamp = IMISSING DO tt = 0,n_times CALL geth_newdate(next_date, trim(start_date), tt*interval_seconds) read(next_date(1:4),*) year read(next_date(6:7),*) month read(next_date(9:10),*) day read(next_date(12:13),*) hours read(next_date(15:16),*) minutes read(next_date(18:19),*) seconds timestamp(tt) = hours*10000 + minutes*100 + seconds if(month >= 2) day = day+31 ! add january if(month >= 3) day = day+28 ! add february if(month >= 4) day = day+31 ! add march if(month >= 5) day = day+30 ! add april if(month >= 6) day = day+31 ! add may if(month >= 7) day = day+30 ! add june if(month >= 8) day = day+31 ! add july if(month >= 9) day = day+31 ! add august if(month >= 10) day = day+30 ! add september if(month >= 11) day = day+31 ! add october if(month >= 12) day = day+30 ! add november if((month > 2) .and. (mod(year,4) == 0)) day = day+1 ! get leap year day datestamp(tt) = (year)*1000 + day END DO END SUBROUTINE v5d_times !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: v5d_proj ! Purpose: Get map information needed to create v5d output file !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE v5d_proj (date, projection, proj_args) USE input_module IMPLICIT NONE character (len=19) :: date integer :: projection real, dimension(100) :: proj_args real :: xlat_s, xlat_n, xlon_w, xlon_e, xlat1, xlon1 real :: xx, yy, temp integer :: istatus, wrftype real, pointer, dimension(:,:,:) :: real_array character (len=3) :: memorder character (len=128) :: cname, stagger character (len=128), dimension(3) :: dimnames integer, dimension(3) :: domain_start, domain_end type (proj_info) :: projst real :: msf real :: hemiflag real :: pi, degtorad domain_start = 1 domain_end(1) = west_east_dim domain_end(2) = south_north_dim domain_end(3) = 1 memorder = "XY" stagger = '-' wrftype = 104 istatus = 0 IF ( iprogram .ge. 6 ) THEN CALL read_spec_field(domain_start, domain_end, 'XLAT', wrftype, & memorder, stagger, dimnames, real_array, date, istatus) ELSE CALL read_spec_field(domain_start, domain_end, 'XLAT_M', wrftype, & memorder, stagger, dimnames, real_array, date, istatus) END IF IF ( istatus == 0) THEN xlat_s = real_array(1,1,1) xlat_n = real_array(1,south_north_dim,1) xlat1 = real_array(1,1,1) END IF istatus = 0 IF ( iprogram .ge. 6 ) THEN CALL read_spec_field(domain_start, domain_end, 'XLONG', wrftype, & memorder, stagger, dimnames, real_array, date, istatus) ELSE CALL read_spec_field(domain_start, domain_end, 'XLONG_M', wrftype, & memorder, stagger, dimnames, real_array, date, istatus) END IF IF ( istatus == 0) THEN xlon_w = real_array(1,south_north_dim,1) xlon_e = real_array(west_east_dim,south_north_dim,1) xlon1 = real_array(1,1,1) END IF IF ( MAP_PROJ == 0 ) THEN !!! ideal cases projection = 0 proj_args(1) = 1. proj_args(2) = 1. proj_args(3) = 1. proj_args(4) = 1. ELSE IF ( MAP_PROJ == 1 ) then !!! Lambert projection = 2 IF ( TRUELAT1 < TRUELAT2 ) THEN !!! truelat1 must be bigger than truelat2 temp = TRUELAT1 TRUELAT1 = TRUELAT2 TRUELAT2 = temp END IF CALL make_lambert(xlat_n,xlon_w,xx,yy) proj_args(1) = TRUELAT1 proj_args(2) = TRUELAT2 proj_args(3) = -yy proj_args(4) = -xx proj_args(5) = -STAND_LON proj_args(6) = 0.001*DX ELSE IF ( MAP_PROJ == 2 ) THEN !!! PS CALL map_set(PROJ_PS,xlat1,xlon1,DX,STAND_LON,TRUELAT1,TRUELAT2, & west_east_dim, south_north_dim, projst) projection = 3 IF (TRUELAT1 .ge. 0.) THEN proj_args(1) = 90. hemiflag = 1. ELSE proj_args(1) = -90. hemiflag = -1. END IF proj_args(2) = -STAND_LON proj_args(3) = float(south_north_dim) - projst%polej proj_args(4) = projst%polei msf = (1.+hemiflag*sin(projst%TRUELAT1*degtorad))/ & (1.+hemiflag*sin(proj_args(1)*degtorad)) proj_args(5) = DX* 0.001/msf ELSE IF ( MAP_PROJ == 3 ) THEN !!! Mercator projection = 1 proj_args(1) = xlat_n proj_args(2) = -xlon_w proj_args(3) = (xlon_w - xlon_e) / (west_east_dim - 1) proj_args(4) = (xlat_n - xlat_s) / (south_north_dim - 1) END IF END SUBROUTINE v5d_proj SUBROUTINE make_lambert(xlatn,xlonw,xx,yy) ! imported from TOVIS5D for MM5 USE input_module real :: xlatn, xlonw, xx, yy real :: radius, degtorad, alpha real :: theta, theta1, theta2, htheta, htheta1, htheta2, cone, dis, rad data radius /6371./ ! put map parameters in as defined by vis5d; ! south is positive direction, so negative ypole pi = 4.0*atan( 1.0 ) degtorad = pi/180.0 dis = 0.001*DX theta = (90.0-xlatn) * degtorad theta1 = (90.0-TRUELAT1) * degtorad theta2 = (90.0-TRUELAT2) * degtorad htheta = theta /2.0 htheta1 = theta1/2.0 htheta2 = theta2/2.0 cone = alog(sin(theta1)/sin(theta2))/alog(tan(htheta1)/tan(htheta2)) alpha = cone * (xlonw - STAND_LON) * degtorad rad = radius * sin(theta1) * (tan(htheta)/tan(htheta1))**cone/cone xx = rad * sin(alpha) / dis yy = rad * cos(alpha) / dis END SUBROUTINE make_lambert END MODULE v5d_module