program get_stddev ! ********************************************************************** ! * . . . ! * PROGRAM: get_stddev ! * ! * PRGMMR: M. ZUPANSKI, D. ZUPANSKI ORG: CIRA/CSU DATE: 2003-09-12 ! * ! * ABSTRACT: THIS SUBROUTINE ESTIMATES THE STANDARD DEVIATION (squre-root variance) ! * AND CONVERTS IT FROM N-DIMENSIONAL VARIABLE TO 1-DIMENSIONAL VARIABLE ! * ! * == Currently, no spatial variability ! * == (Just assign the standard deviation value for each variable) ! * ! * INPUT: ! * cvar%stddev .... parameter defining an estimate of the standard deviation ! * bias_factor .... defines reduction of magnitude of bias with respect to ! * the initial conditions standard deviation ! * ! * OUTPUT: ! * stddev .... standard deviation to be used in pertubing random perturbations ! * in cold start of ensemble data assimilation algorithm ! * ! * PROGRAM LOG: ! * ! * 08/14/2003 ..... M. ZUPANSKI, D. ZUPANSKI ! * 09/12/2003 ..... D. ZUPANSKI ! * 09/12/2003 ..... M. ZUPANSKI: Add standard deviation estimate to the 'convert' program ! * ! ------------------------------------------------------------------------- ! * Keep these functions for more complex options in future: ! * == Temperature, Wind, ...: y=alf*exp[-bet*((p-pref)/1000)**2] ! * == Moisture, Mixing ratio: y=alf*(1-exp[-bet*((p-pref)/1000)**2]) ! ********************************************************************** real,parameter :: bias_factor=5.0E-02 ! real,parameter :: bias_factor=1.0E-04 ! real,parameter :: bias_factor=1.0E-02 integer,parameter :: imodl=20 ! model dimension namelist integer,parameter :: icntrl=21 ! input control variable file integer,parameter :: ix1d=22 ! 1-d control variable integer,parameter :: idev_input=25 ! input standard deviation integer,parameter :: idev=105 ! output standard deviation INCLUDE 'define_var_list_type.h' INCLUDE 'define_cvar_type.h' TYPE (var_list_type), dimension(:), allocatable :: cvar_list TYPE (cvar_type), dimension(:), allocatable :: cvar real, dimension(:), allocatable :: stddev,stddev_input integer, dimension(1) :: NNXP,NNYP,NNZP,NNSOIL,NNSNOW integer :: cvar_max integer :: nbias,num_bias character*4 :: chmember character*20 :: filename character*9 :: cv_name,cv_name_read logical :: open_ic,open_param,open_bias !------------------------------ NAMELIST /MODEL_DIMENSION/ NNXP,NNYP,NNZP,NNSOIL,NNSNOW !------------------------------ write(*,*) " start get_stddev" !--- get model dimensions write(filename,200) 200 format('model_dimension.name') CLOSE(imodl) OPEN(UNIT=imodl,FILE=filename,FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(3,*) imodl,' OPEN UNIT ERROR IER=',IER READ(imodl,MODEL_DIMENSION) CLOSE(imodl,STATUS='KEEP') write(*,*) "NNXP=",NNXP(1)," NNYP=",NNYP(1)," NNZP=",NNZP(1) INCLUDE 'cntrl_vrbl_list.h' write(*,*) "max_num_of_cntrl_vrbls=",max_num_of_cntrl_vrbls do ncv=1,max_num_of_cntrl_vrbls write(*,*) "ndim ",cvar_list(ncv)%ndim write(*,*) "start1 ",cvar_list(ncv)%start_index(1) write(*,*) "start2 ",cvar_list(ncv)%start_index(2) write(*,*) "start3 ",cvar_list(ncv)%start_index(3) write(*,*) "start4 ",cvar_list(ncv)%start_index(4) write(*,*) "end1 ",cvar_list(ncv)%end_index(1) write(*,*) "end2 ",cvar_list(ncv)%end_index(2) write(*,*) "end3 ",cvar_list(ncv)%end_index(3) write(*,*) "end4 ",cvar_list(ncv)%end_index(4) write(*,*) "name ",cvar_list(ncv)%name write(*,*) "std dev ",cvar_list(ncv)%stddev write(*,*) "description ",cvar_list(ncv)%description end do CLOSE(icntrl) OPEN(icntrl,FILE='cntrl_vrbl',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) icntrl,' OPEN UNIT ERROR IER=',IER READ(icntrl,*) cvar_max allocate(cvar(1:cvar_max)) do i=1,cvar_max READ(icntrl,*) cvar(i) write(*,*) cvar(i) if(cvar(i)%bias) then num_bias=cvar(i)%num_bias else num_bias=0 end if end do CLOSE(icntrl,STATUS='KEEP') !========== start reading 1d control variable ============== ! open input files write(filename,1000) 1000 format('x1d_cntrl') write(*,*) "input: 1d filename=",filename CLOSE(ix1d) OPEN(UNIT=ix1d,FILE=filename,FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) ix1d,' OPEN UNIT ERROR IER=',IER READ(ix1d) idim_1d READ(ix1d) CLOSE(ix1d,status='KEEP') write(filename,1020) 1020 format('fcst_std_dev_input') write(*,*) "input: fcst_std_dev_input filename=",filename CLOSE(idev_input) OPEN(UNIT=idev_input,FILE=filename,FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) idev_input,' OPEN UNIT ERROR IER=',IER allocate(stddev(1:idim_1d)) !---------------------------------------------- ! 100 format(E20.10) icount=0 nvar=0 do i=1,cvar_max write(*,*) "do i=",i do j=1,max_num_of_cntrl_vrbls write(*,*) "do j=",j if(cvar(i)%name.eq.cvar_list(j)%name) then nvar=nvar+1 i1s=cvar_list(j)%start_index(1) i2s=cvar_list(j)%start_index(2) i3s=cvar_list(j)%start_index(3) i4s=cvar_list(j)%start_index(4) i1e=cvar_list(j)%end_index(1) i2e=cvar_list(j)%end_index(2) i3e=cvar_list(j)%end_index(3) i4e=cvar_list(j)%end_index(4) !-- initial conditions if(cvar(i)%ic) then if(cvar_list(j)%ndim.eq.4) then do i4=i4s,i4e do i3=i3s,i3e do i2=i2s,i2e do i1=i1s,i1e icount=icount+1 stddev(icount)=cvar_list(j)%stddev end do end do end do end do else if(cvar_list(j)%ndim.eq.3) then allocate(stddev_input(i3s:i3e)) rewind idev_input do read(idev_input,*,end=220) cv_name read(idev_input,*) stddev_input if(cv_name.eq.cvar(i)%name) then exit else cycle end if enddo if(cv_name.eq.'u ') then stddev_input(:)=cvar_list(j)%stddev ! stddev_input(:)=stddev_input(:)*1.0 ! stddev_input(:)=stddev_input(:)*0.5 else if(cv_name.eq.'v ') then stddev_input(:)=cvar_list(j)%stddev ! stddev_input(:)=stddev_input(:)*1.0 ! stddev_input(:)=stddev_input(:)*0.5 else if(cv_name.eq.'w ') then stddev_input(:)=cvar_list(j)%stddev ! stddev_input(:)=stddev_input(:)*1.0 ! stddev_input(:)=stddev_input(:)*0.5 else if(cv_name.eq.'exnr ') then stddev_input(:)=cvar_list(j)%stddev ! stddev_input(:)=stddev_input(:)*1.0 ! stddev_input(:)=stddev_input(:)*0.5 else if(cv_name.eq.'thetail ') then stddev_input(:)=cvar_list(j)%stddev ! stddev_input(:)=stddev_input(:)*1.0 ! stddev_input(:)=stddev_input(:)*0.5 else if(cv_name.eq.'r_total ') then ! stddev_input(:)=cvar_list(j)%stddev stddev_input(:)=stddev_input(:)*1.0 ! stddev_input(:)=stddev_input(:)*0.5 endif write(*,*) " cv_name=",cv_name write(*,*) " stddev_input=",stddev_input do i3=i3s,i3e do i2=i2s,i2e do i1=i1s,i1e icount=icount+1 ! stddev(icount)=cvar_list(j)%stddev stddev(icount)=stddev_input(i3) end do end do end do deallocate(stddev_input) else if(cvar_list(j)%ndim.eq.2) then do i2=i2s,i2e do i1=i1s,i1e icount=icount+1 stddev(icount)=cvar_list(j)%stddev end do end do else if(cvar_list(j)%ndim.eq.1) then do i1=i1s,i1e icount=icount+1 stddev(icount)=cvar_list(j)%stddev end do else write(*,*) "IC DOES NOT HAVE APPROPRIATE ndim=",cvar_list(j)%ndim end if end if !-- empirical model parameters if(cvar(i)%param) then icount=icount+1 stddev(icount)=cvar_list(j)%stddev end if !-- model bias if(cvar(i)%bias) then do nbias=1,num_bias if(cvar_list(j)%ndim.eq.4) then do i4=i4s,i4e do i3=i3s,i3e do i2=i2s,i2e do i1=i1s,i1e icount=icount+1 stddev(icount)=bias_factor*cvar_list(j)%stddev end do end do end do end do else if(cvar_list(j)%ndim.eq.3) then do i3=i3s,i3e do i2=i2s,i2e do i1=i1s,i1e icount=icount+1 stddev(icount)=bias_factor*cvar_list(j)%stddev end do end do end do else if(cvar_list(j)%ndim.eq.2) then do i2=i2s,i2e do i1=i1s,i1e icount=icount+1 stddev(icount)=bias_factor*cvar_list(j)%stddev end do end do else if(cvar_list(j)%ndim.eq.1) then do i1=i1s,i1e icount=icount+1 stddev(icount)=bias_factor*cvar_list(j)%stddev end do else write(*,*) "BIAS DOES NOT HAVE APPROPRIATE ndim=",cvar_list(j)%ndim end if end do end if !---------------------------------------------- end if end do end do !---------------------------------------------- if(icount.ne.idim_1d) then write(*,*) " get_stddev, PROBLEM with icount=",icount, & " idim_1d=",idim_1d stop endif CLOSE(idev_input,status='KEEP') !---------write 1d standard deviation file------------------------- !-- read standard deviation CLOSE(idev) OPEN(idev,file='std_dev',status='unknown',form='unformatted') WRITE(idev) idim_1d WRITE(idev) stddev CLOSE(idev,status='KEEP') write(*,*) " end get_stddev" stop 220 write(*,*) " PROBLEM: cv_name on idev_input is not appropriate" write(*,*) " end get_stddev" stop end program get_stddev