!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !====================================================================== ! BEGINNING OF ENSDA_INIT !====================================================================== subroutine ensda_init & (n,dt_model,my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi & ,visc_model) ! ********************************************************************** ! ! PROGRAM: ensda_init read and write initial ensemble data assimilation fields ! PRGMMR: M. ZUPANSKI ORG: CIRA/CSU DATE: 2004-05-10 ! ! ! REVISION HISTORY: ! ! 05/10/2004 ..... M. ZUPANSKI ! !----------------------------------------------------------------------- use ensda_variables implicit NONE real(kind=selected_real_kind(13)) :: dt_model real(kind=selected_real_kind(13)) :: visc_model integer :: n integer :: my_task integer :: im,jm,iim,jjm,nsdm real(kind=selected_real_kind(13)),dimension(iim,jjm,nsdm) :: h,psi,chi character(*) :: swm_name real :: visc !------- visc = parameter read in namelists !------- visc_model = parameter used by model (double precision in this case) !------- dt_model = parameter used by model (double precision in this case) !=========================================================== nbias=0 num_bias=0 open_ic=.FALSE. open_param=.FALSE. open_bias=.FALSE. call read_write_ensda_namelists & (my_task,im,jm,iim,jjm,nsdm,visc) visc_model=visc !----------- time=int(fcst_length/(dt_model/60.)) !----------- if(my_task.eq.0) then write(*,*) "ensda_init: n=",n write(*,*) "ensda_init: dt_model=",dt_model write(*,*) "ensda_init: swm_name=",swm_name write(*,*) "ensda_init: im,jm,iim,jjm,nsdm=",im,jm,iim,jjm,nsdm write(*,*) "ensda_init: h min,max=",minval(h),maxval(h) write(*,*) "ensda_init: psi min,max=",minval(psi),maxval(psi) write(*,*) "ensda_init: chi min,max=",minval(chi),maxval(chi) write(*,*) "ensda_init: visc=",visc write(*,*) "ensda_init: visc_model=",visc_model end if call write_history_format & (n,dt_model,my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi) call write_history_nonformat & (my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi) call write_model_ic & (my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi) call write_model_param & (my_task,visc) do while(nbias.le.num_bias) call write_model_bias & (my_task,swm_name,im,jm,iim,jjm,nsdm) nbias=nbias+1 enddo end subroutine ensda_init !====================================================================== ! END OF ENSDA_INIT !====================================================================== !====================================================================== ! BEGINNING OF INITIALIZE_ENSDA !====================================================================== subroutine initialize_ensda & (n,dt_model,my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi & ,visc_model) ! ********************************************************************** ! ! PROGRAM: ensda_init read and initialize ensemble data assimilation fields ! PRGMMR: M. ZUPANSKI ORG: CIRA/CSU DATE: 2004-05-10 ! ! ! REVISION HISTORY: ! ! 05/10/2004 ..... M. ZUPANSKI ! !----------------------------------------------------------------------- use ensda_variables implicit NONE real(kind=selected_real_kind(13)) :: dt_model real(kind=selected_real_kind(13)) :: visc_model integer :: n integer :: my_task integer :: im,jm,iim,jjm,nsdm real(kind=selected_real_kind(13)),dimension(iim,jjm,nsdm) :: h,psi,chi character(*) :: swm_name real :: visc !------- visc = parameter read in namelists !------- visc_model = parameter used by model (double precision in this case) !------- dt_model = parameter used by model (double precision in this case) !=========================================================== nbias=0 num_bias=0 open_ic=.FALSE. open_param=.FALSE. open_bias=.FALSE. call read_ensda_namelists & (my_task,im,jm,iim,jjm,nsdm,visc) !----------- time=int(fcst_length/(dt_model/60.)) !----------- call read_history_nonformat & (my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi) call read_model_ic & (my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi) call read_model_param & (my_task,visc) !-- assign model viscosity from input namelist --------- visc_model=visc if(my_task.eq.0) then write(*,*) "initialize_ensda: n=",n write(*,*) "initialize_ensda: dt_model=",dt_model write(*,*) "initialize_ensda: time=",time write(*,*) "initialize_ensda: fcst_length=",fcst_length write(*,*) "initialize_ensda: swm_name=",swm_name write(*,*) "initialize_ensda: im,jm,iim,jjm,nsdm=",im,jm,iim,jjm,nsdm write(*,*) "initialize_ensda: h min,max=",minval(h),maxval(h) write(*,*) "initialize_ensda: psi min,max=",minval(psi),maxval(psi) write(*,*) "initialize_ensda: chi min,max=",minval(chi),maxval(chi) write(*,*) "initialize_ensda: visc=",visc write(*,*) "initialize_ensda: visc_model=",visc_model end if !-- read _phi as _bias (need to rename it after read) call read_model_bias & (my_task,swm_name,im,jm,iim,jjm,nsdm) if(num_bias.gt.0) then nerr=time/max(1,num_bias) else nerr=99999 end if if(associated(h_phi_ptr)) then h_phi(:,:,:)=h_bias(:,:,:) h_bias(:,:,:)=0.0 end if if(associated(psi_phi_ptr)) then psi_phi(:,:,:)=psi_bias(:,:,:) psi_bias(:,:,:)=0.0 end if if(associated(chi_phi_ptr)) then chi_phi(:,:,:)=chi_bias(:,:,:) chi_bias(:,:,:)=0.0 end if end subroutine initialize_ensda !====================================================================== ! END OF INITIALIZE_ENSDA !====================================================================== !====================================================================== ! BEGINNING OF READ_ENSDA_NAMELISTS !====================================================================== subroutine read_ensda_namelists & (my_task,im,jm,iim,jjm,nsdm,visc) ! ********************************************************************** ! * PROGRAM: read_ensda_namelists ! * ! * PRGMMR: D. ZUPANSKI ORG: CIRA/CSU DATE: 2003-12-12 ! * ! * ABSTRACT: Read namelists and other constant input files for use in ensda ! ! ! * PROGRAM LOG: ! * ! * 12/12/2003 ..... D. ZUPANSKI, M. ZUPANSKI ! * 05/06/2004 ..... M. ZUPANSKI: Adjust to SWM_CSU model use !==================================== use ensda_variables implicit none integer :: my_task integer :: im,jm,iim,jjm,nsdm integer :: imax_gbl,jmax_gbl real :: visc integer :: IER,i integer,dimension(:),allocatable :: & cvar_num_bias ! number of biases per cntrl vrbl character (len=9),dimension(:),allocatable :: & cvar_name ! cntrl vrbl name logical,dimension(:),allocatable :: & cvar_ic, &! init cond logical mask (per vrbl) cvar_param, &! empir param logical mask (per vrbl) cvar_bias ! model bias logical mask (per vrbl) !------------------------------ NAMELIST /MODEL_DIMENSION/ NNXP,NNYP,NNZP,NNSOIL,NNSNOW NAMELIST /LENGTH/ fcst_length NAMELIST /MODEL_CNSTS/ visc,aa !===================================== !--- get model dimensions CLOSE(500) OPEN(UNIT=500,FILE='model_dimension.name',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*)' 500 OPEN UNIT ERROR IER=',IER READ(500,MODEL_DIMENSION) CLOSE(500,STATUS='KEEP') if(my_task.eq.0) then write(*,*) "NNXP=",NNXP(1)," NNYP=",NNYP(1)," NNZP=",NNZP(1) end if imax_gbl=NNXP(1) jmax_gbl=NNYP(1) if(my_task.eq.0) then if(im.eq.imax_gbl.and.jm.eq.jmax_gbl) then write(*,*) "GLOBAL DIMENSIONS ARE CORRECT: im,jm=",im,jm else write(*,*) "GLOBAL DIMENSIONS ARE NOT CORRECT: im,jm=",im,jm write(*,*) " imax_gbl,jmax_gbl=",imax_gbl,jmax_gbl stop end if end if !----------------------------------------- !----read cnsts file-------------------------------------------------- CLOSE(501) OPEN(UNIT=501,FILE='model_cnsts.name',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*)' 501 OPEN UNIT ERROR IER=',IER read(501,MODEL_CNSTS) if(my_task.eq.0) then if(visc.lt.0.) then write(*,*) "PROBLEM!!! NEGATIVE VISCOUSITY PARAMETER visc=",visc write(*,*) "CHANGE VISCOSITY PARAMETER TO visc=0.0" visc=0.0 else write(*,*) "CURRENT VISCOSITY PARAMETER visc=",visc end if end if !------------------------------------------ CLOSE(502) OPEN(UNIT=502,FILE='cycle_date',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*) ' 502 OPEN UNIT ERROR IER=',IER read(502,300) cycle_date 300 format(a14) !--- get fcst_length CLOSE(503) OPEN(UNIT=503,FILE='fcst_length.name',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*)' 503 OPEN UNIT ERROR IER=',IER READ(503,LENGTH) CLOSE(503,STATUS='KEEP') if(my_task.eq.0) then write(*,*) "fcst_length=",fcst_length," minutes" end if !--- read control variable list ------- CLOSE(103) OPEN(UNIT=103,FILE='cntrl_vrbl',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*)' 103 OPEN UNIT ERROR IER=',IER READ(103,*) cvar_max !-- read control variable parameters --------------- allocate(cvar_name(1:cvar_max)) allocate(cvar_ic(1:cvar_max)) allocate(cvar_param(1:cvar_max)) allocate(cvar_bias(1:cvar_max)) allocate(cvar_num_bias(1:cvar_max)) nullify(h_phi_ptr) nullify(h_bias_ptr) nullify(psi_phi_ptr) nullify(psi_bias_ptr) nullify(chi_phi_ptr) nullify(chi_bias_ptr) do i=1,cvar_max READ(103,*) cvar_name(i),cvar_ic(i),cvar_param(i),cvar_bias(i),cvar_num_bias(i) open_ic=open_ic.or.cvar_ic(i) open_param=open_param.or.cvar_param(i) open_bias=open_bias.or.cvar_bias(i) if(cvar_bias(i)) then if(cvar_name(i).eq.'h ') then allocate(h_phi(iim,jjm,nsdm)) allocate(h_bias(iim,jjm,nsdm)) h_phi(:,:,:)=0.0 h_bias(:,:,:)=0.0 h_phi_ptr => h_phi h_bias_ptr => h_bias else if(cvar_name(i).eq.'psi ') then allocate(psi_phi(iim,jjm,nsdm)) allocate(psi_bias(iim,jjm,nsdm)) psi_phi(:,:,:)=0.0 psi_bias(:,:,:)=0.0 psi_phi_ptr => psi_phi psi_bias_ptr => psi_bias else if(cvar_name(i).eq.'chi ') then allocate(chi_phi(iim,jjm,nsdm)) allocate(chi_bias(iim,jjm,nsdm)) chi_phi(:,:,:)=0.0 chi_bias(:,:,:)=0.0 chi_phi_ptr => chi_phi chi_bias_ptr => chi_bias else write(*,*)"WARNING: NO model bias found for i=",i end if !!!! if(cvar_name(i).eq.'h ') end if !!!! if(cvar_bias(i)) then end do !!!! do i=1,cvar_max if(open_bias) then num_bias=maxval(cvar_num_bias) else num_bias=0 end if deallocate(cvar_name) deallocate(cvar_ic) deallocate(cvar_param) deallocate(cvar_bias) deallocate(cvar_num_bias) CLOSE(103,STATUS='KEEP') !------------------------------------------------- end subroutine read_ensda_namelists !====================================================================== ! END OF READ_ENSDA_NAMELISTS !====================================================================== !====================================================================== ! BEGINNING OF READ_WRITE_ENSDA_NAMELISTS !====================================================================== subroutine read_write_ensda_namelists & (my_task,im,jm,iim,jjm,nsdm,visc) ! ********************************************************************** ! * PROGRAM: write_ensda_namelists ! * ! * PRGMMR: D. ZUPANSKI ORG: CIRA/CSU DATE: 2003-12-12 ! * ! * ABSTRACT: Read namelists and other constant input files for use in ensda ! ! ! * PROGRAM LOG: ! * ! * 12/12/2003 ..... D. ZUPANSKI, M. ZUPANSKI ! * 05/06/2004 ..... M. ZUPANSKI: Adjust to SWM_CSU model use !==================================== use ensda_variables implicit none integer :: my_task integer :: im,jm,iim,jjm,nsdm integer :: imax_gbl,jmax_gbl real :: visc integer :: IER,i integer,dimension(:),allocatable :: & cvar_num_bias ! number of biases per cntrl vrbl character (len=9),dimension(:),allocatable :: & cvar_name ! cntrl vrbl name logical,dimension(:),allocatable :: & cvar_ic, &! init cond logical mask (per vrbl) cvar_param, &! empir param logical mask (per vrbl) cvar_bias ! model bias logical mask (per vrbl) !------------------------------ NAMELIST /MODEL_CNSTS/ visc,aa !===================================== !-- read cycle_date ----------------------- CLOSE(502) OPEN(UNIT=502,FILE='cycle_date',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*) ' 502 OPEN UNIT ERROR IER=',IER read(502,300) cycle_date 300 format(a14) !--- write model dimensions NNXP(1)=im NNYP(1)=jm NNZP(1)=1 NNSOIL(1)=1 NNSNOW(1)=1 if(my_task.eq.0) then CLOSE(105) OPEN(UNIT=105,FILE='model_dimension',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*)' 105 OPEN UNIT ERROR IER=',IER WRITE(105,1001) NNXP(:) WRITE(105,1002) NNYP(:) WRITE(105,1003) NNZP(:) WRITE(105,1004) NNSOIL(:) WRITE(105,1005) NNSNOW(:) CLOSE(105,STATUS='KEEP') 1001 format(' NNXP = ',I8) 1002 format(' NNYP = ',I8) 1003 format(' NNZP = ',I8) 1004 format(' NNSOIL = ',I8) 1005 format(' NNSNOW = ',I8) !-- if nested, number of nests=n, it would assume NNXP(1:n) ---------- !1001 format(' NNXP = ',nI8) !1002 format(' NNYP = ',nI8) !1003 format(' NNZP = ',nI8) !1004 format(' NNSOIL = ',nI8) !1005 format(' NNSNOW = ',nI8) !-------------------------------------------------------------------- write(*,*) "NNXP=",NNXP(:)," NNYP=",NNYP(:)," NNZP=",NNZP(:) end if !----------------------------------------- !----read cnsts file-------------------------------------------------- CLOSE(501) OPEN(UNIT=501,FILE='model_cnsts.name',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*)' 501 OPEN UNIT ERROR IER=',IER read(501,MODEL_CNSTS) if(my_task.eq.0) then if(visc.lt.0.) then write(*,*) "PROBLEM!!! NEGATIVE VISCOUSITY PARAMETER visc=",visc write(*,*) "CHANGE VISCOSITY PARAMETER TO visc=0.0" visc=0.0 else write(*,*) "CURRENT VISCOSITY PARAMETER visc=",visc end if end if !------------------------------------------ !--- read control variable list ------- CLOSE(103) OPEN(UNIT=103,FILE='cntrl_vrbl',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*)' 103 OPEN UNIT ERROR IER=',IER READ(103,*) cvar_max !-- read control variable parameters --------------- allocate(cvar_name(1:cvar_max)) allocate(cvar_ic(1:cvar_max)) allocate(cvar_param(1:cvar_max)) allocate(cvar_bias(1:cvar_max)) allocate(cvar_num_bias(1:cvar_max)) nullify(h_phi_ptr) nullify(h_bias_ptr) nullify(psi_phi_ptr) nullify(psi_bias_ptr) nullify(chi_phi_ptr) nullify(chi_bias_ptr) do i=1,cvar_max READ(103,*) cvar_name(i),cvar_ic(i),cvar_param(i),cvar_bias(i),cvar_num_bias(i) open_ic=open_ic.or.cvar_ic(i) open_param=open_param.or.cvar_param(i) open_bias=open_bias.or.cvar_bias(i) if(cvar_bias(i)) then if(cvar_name(i).eq.'h ') then allocate(h_phi(iim,jjm,nsdm)) allocate(h_bias(iim,jjm,nsdm)) h_phi(:,:,:)=0.0 h_bias(:,:,:)=0.0 h_phi_ptr => h_phi h_bias_ptr => h_bias else if(cvar_name(i).eq.'psi ') then allocate(psi_phi(iim,jjm,nsdm)) allocate(psi_bias(iim,jjm,nsdm)) psi_phi(:,:,:)=0.0 psi_bias(:,:,:)=0.0 psi_phi_ptr => psi_phi psi_bias_ptr => psi_bias else if(cvar_name(i).eq.'chi ') then allocate(chi_phi(iim,jjm,nsdm)) allocate(chi_bias(iim,jjm,nsdm)) chi_phi(:,:,:)=0.0 chi_bias(:,:,:)=0.0 chi_phi_ptr => chi_phi chi_bias_ptr => chi_bias else write(*,*)"WARNING: NO model bias found for i=",i end if !!!! if(cvar_name(i).eq.'h ') end if !!!! if(cvar_bias(i)) then end do !!!! do i=1,cvar_max if(open_bias) then num_bias=maxval(cvar_num_bias) else num_bias=0 end if deallocate(cvar_name) deallocate(cvar_ic) deallocate(cvar_param) deallocate(cvar_bias) deallocate(cvar_num_bias) CLOSE(103,STATUS='KEEP') !------------------------------------------------- end subroutine read_write_ensda_namelists !====================================================================== ! END OF READ_WRITE_ENSDA_NAMELISTS !====================================================================== !====================================================================== ! BEGINNING OF READ_HISTORY_NONFORMAT !====================================================================== subroutine read_history_nonformat & (my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi) ! ********************************************************************** ! * PROGRAM: read_history ! * ! * PRGMMR: D. ZUPANSKI ORG: CIRA/CSU DATE: 2003-12-12 ! * ! * ABSTRACT: Read history file (non-formatted) ! ! ! * PROGRAM LOG: ! * ! * 12/12/2003 ..... D. ZUPANSKI, M. ZUPANSKI ! * 05/06/2004 ..... M. ZUPANSKI: Adjust to SWM_CSU model use !==================================== use gather_scatter use ensda_variables implicit none integer :: my_task integer :: im,jm,iim,jjm,nsdm real (kind=selected_real_kind(13)),dimension(iim,jjm,nsdm) :: h,psi,chi character(*) :: swm_name integer :: IER !------------------------------------- !----read history file----------------------------------------------- allocate(h_gbl(1:im,1:jm)) allocate(psi_gbl(1:im,1:jm)) allocate(chi_gbl(1:im,1:jm)) CLOSE(400) OPEN(UNIT=400,FILE='model_history',FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*)' 400 OPEN UNIT ERROR IER=',IER read(400) h_gbl read(400) psi_gbl read(400) chi_gbl CLOSE(400,STATUS='KEEP') allocate(ensda_work_gbl(1:im,1:jm)) ensda_work_gbl(:,:)=h_gbl(:,:) call scatter(swm_name,im,jm,iim,jjm,1,nsdm,ensda_work_gbl,h) deallocate(h_gbl) ensda_work_gbl(:,:)=psi_gbl(:,:) call scatter(swm_name,im,jm,iim,jjm,1,nsdm,ensda_work_gbl,psi) deallocate(psi_gbl) ensda_work_gbl(:,:)=chi_gbl(:,:) call scatter(swm_name,im,jm,iim,jjm,1,nsdm,ensda_work_gbl,chi) deallocate(chi_gbl) deallocate(ensda_work_gbl) !------------------------------------- end subroutine read_history_nonformat !====================================================================== ! END OF READ_HISTORY_NONFORMAT !====================================================================== !====================================================================== ! BEGINNING OF WRITE_HISTORY_NONFORMAT !====================================================================== subroutine write_history_nonformat & (my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi) ! ********************************************************************** ! * PROGRAM: write_history ! * ! * PRGMMR: D. ZUPANSKI ORG: CIRA/CSU DATE: 2003-12-12 ! * ! * ABSTRACT: Read history file (non-formatted) ! ! ! * PROGRAM LOG: ! * ! * 12/12/2003 ..... D. ZUPANSKI, M. ZUPANSKI ! * 05/06/2004 ..... M. ZUPANSKI: Adjust to SWM_CSU model use !==================================== use gather_scatter use ensda_variables implicit none integer :: my_task integer :: im,jm,iim,jjm,nsdm real (kind=selected_real_kind(13)),dimension(iim,jjm,nsdm) :: h,psi,chi character(*) :: swm_name integer :: IER !------------------------------------- allocate(ensda_work_gbl(1:im,1:jm)) call gather(swm_name,im,jm,iim,jjm,1,nsdm,h,ensda_work_gbl) allocate(h_gbl(1:im,1:jm)) h_gbl(:,:)=ensda_work_gbl(:,:) call gather(swm_name,im,jm,iim,jjm,1,nsdm,psi,ensda_work_gbl) allocate(psi_gbl(1:im,1:jm)) psi_gbl(:,:)=ensda_work_gbl(:,:) call gather(swm_name,im,jm,iim,jjm,1,nsdm,chi,ensda_work_gbl) allocate(chi_gbl(1:im,1:jm)) chi_gbl(:,:)=ensda_work_gbl(:,:) deallocate(ensda_work_gbl) !----write history file----------------------------------------------- if(my_task.eq.0) then CLOSE(400) OPEN(UNIT=400,FILE='model_history',FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*)' 400 OPEN UNIT ERROR IER=',IER write(400) h_gbl write(400) psi_gbl write(400) chi_gbl CLOSE(400,STATUS='KEEP') end if deallocate(h_gbl) deallocate(psi_gbl) deallocate(chi_gbl) !------------------------------------- end subroutine write_history_nonformat !====================================================================== ! END OF WRITE_HISTORY_NONFORMAT !====================================================================== !====================================================================== ! BEGINNING OF WRITE_HISTORY_FORMAT !====================================================================== subroutine write_history_format & (n,dt_sec,my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi) ! ********************************************************************** ! * PROGRAM: write_history_format ! * ! * PRGMMR: D. ZUPANSKI ORG: CIRA/CSU DATE: 2003-12-12 ! * ! * ABSTRACT: Write history file (formatted) ! ! ! * PROGRAM LOG: ! * ! * 12/12/2003 ..... D. ZUPANSKI, M. ZUPANSKI ! * 05/06/2004 ..... M. ZUPANSKI: Adjust to SWM_CSU model use !==================================== use gather_scatter use ensda_variables implicit none integer :: n real(selected_real_kind(13)) :: dt_sec integer :: my_task integer :: im,jm,iim,jjm,nsdm real (kind=selected_real_kind(13)),dimension(iim,jjm,nsdm) :: h,psi,chi character(*) :: swm_name integer,parameter :: iparm=300 ! PARM file unit integer,parameter :: ihist=401 ! formatted history file integer :: IER,i,j character(len=40) :: filename character(len=14) :: rdate_out character(len=14) :: cycle_start_date integer :: cycle_interval,N_cycles character(len=13) :: jdate,jdate_fg integer :: year,jday,hour integer :: year_fg,jday_fg,hour_fg integer :: idiff,ndiff character(len=9) :: h_name,psi_name,chi_name ! ! DECLARE NAMELISTS ! NAMELIST /CYCLE/ cycle_start_date,cycle_interval,N_cycles !=========================================================== allocate(ensda_work_gbl(1:im,1:jm)) call gather(swm_name,im,jm,iim,jjm,1,nsdm,h,ensda_work_gbl) allocate(h_gbl(1:im,1:jm)) h_gbl(:,:)=ensda_work_gbl(:,:) call gather(swm_name,im,jm,iim,jjm,1,nsdm,psi,ensda_work_gbl) allocate(psi_gbl(1:im,1:jm)) psi_gbl(:,:)=ensda_work_gbl(:,:) call gather(swm_name,im,jm,iim,jjm,1,nsdm,chi,ensda_work_gbl) allocate(chi_gbl(1:im,1:jm)) chi_gbl(:,:)=ensda_work_gbl(:,:) deallocate(ensda_work_gbl) !------------------------------------- ! define names write(h_name,200) 200 format('height') write(psi_name,210) 210 format('psi') write(chi_name,220) 220 format('chi') 100 format(E20.10) 110 format(2i7) 300 format(a9) !----write history file at the end of forecast----------------------- if(my_task.eq.0) then CLOSE(iparm) OPEN(UNIT=iparm,FILE='cycle.parm',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*) iparm,' OPEN UNIT ERROR IER=',IER READ(iparm,CYCLE) CLOSE(iparm,status='KEEP') !----- idiff=(sec), cycle_interval=(min) idiff=cycle_interval*60 ndiff=int(n*dt_sec+0.5) if(MOD(ndiff,idiff).eq.0) then call rdate_to_jdate (cycle_date,jdate,year,jday,hour) call add_to_jdate(year,jday,hour,ndiff,year_fg, & jday_fg,hour_fg) call jdate_make_big (year_fg,jday_fg,hour_fg,jdate_fg) call jdate_to_rdate (jdate_fg,rdate_out) write(filename,1001) rdate_out 1001 format('model_history.',a14) !------------------------------------- CLOSE(ihist) OPEN(UNIT=ihist,FILE=filename,FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) ihist,' OPEN UNIT ERROR IER=',IER write(ihist,110) im,jm write(ihist,300) h_name !-- height do i=1,im do j=1,jm write(ihist,100) h_gbl(i,j) end do end do write(ihist,300) psi_name !-- stream function do i=1,im do j=1,jm write(ihist,100) psi_gbl(i,j) end do end do write(ihist,300) chi_name !-- velocity potential do i=1,im do j=1,jm write(ihist,100) chi_gbl(i,j) end do end do CLOSE(ihist,STATUS='KEEP') end if !!! (my_task=0) deallocate(h_gbl) deallocate(psi_gbl) deallocate(chi_gbl) end if end subroutine write_history_format !====================================================================== ! END OF WRITE_HISTORY_FORMAT !====================================================================== !====================================================================== ! BEGINNING OF WRITE_FG !====================================================================== subroutine write_fg & (n,dt_sec,my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi) ! ********************************************************************** ! ! ROUTINE: write_fg: write files with first guess forecast ! used in ensda ! PRGMMR: D. ZUPANSKI ORG: CIRA/CSU DATE: 2003-09-30 ! ! ! REVISION HISTORY: ! ! 09/30/2003 ..... D. ZUPANSKI ! !----------------------------------------------------------------------- use gather_scatter use ensda_variables implicit NONE integer :: n real(selected_real_kind(13)) :: dt_sec integer :: my_task integer :: im,jm,iim,jjm,nsdm real(kind=selected_real_kind(13)),dimension(iim,jjm,nsdm) :: h,psi,chi character(*) :: swm_name integer,parameter :: i_parm=201 ! PARM file unit integer,parameter :: i_obstype=202 ! OBSTYPE file unit integer,parameter :: i_jdate=203 ! jdate_ file unit integer,parameter :: i_fguess=204 ! first guess file unit integer :: fgvar_max,IER,i character(len=40) :: file_fguess,file_obstype,file_jdate character(len=14) :: rdate_out character(len=13) :: jdate,jdate_fg,jdate_obs character(len=9), dimension(:),allocatable :: fgvar_name character (len=6) :: OBSTYPE character (len=3) :: OBSFLAG integer :: kk integer :: year,jday,hour integer :: year_fg,jday_fg,hour_fg integer :: year_obs,jday_obs,hour_obs integer :: iobs,iobsmax,delobs,idiff,idiff_obs integer, dimension(:),allocatable :: kk_start,kk_end logical :: output ! ! DECLARE NAMELISTS ! NAMELIST /OBSPARM/ OBSTYPE,delobs,OBSFLAG !=========================================================== CLOSE(i_parm) OPEN(UNIT=i_parm,FILE='kk_index.parm',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*) i_parm,' OPEN UNIT ERROR IER=',IER read(i_parm,*) iobsmax allocate(kk_start(1:iobsmax)) allocate(kk_end(1:iobsmax)) do i=1,iobsmax read(i_parm,*) kk_start(i),kk_end(i) end do CLOSE(i_parm,status='keep') output=.FALSE. do iobs=1,iobsmax write(file_obstype,1002) iobs 1002 format('OBSTYPE_',i2.2,'.name') CLOSE(i_obstype) OPEN(UNIT=i_obstype,FILE=file_obstype,FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) i_obstype,' OPEN UNIT ERROR IER=',IER read(i_obstype,OBSPARM) CLOSE(i_obstype,status='KEEP') if(OBSFLAG.eq.'YES') then write(file_fguess,1004) OBSTYPE 1004 format('first_guess_',a6) CLOSE(105) OPEN(105,FILE=file_fguess,FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*)' 105 OPEN UNIT ERROR IER=',IER READ(105,*) fgvar_max allocate(fgvar_name(1:fgvar_max)) do i=1,fgvar_max READ(105,*) fgvar_name(i) end do CLOSE(105,STATUS='KEEP') !----- define appropriate jdate call rdate_to_jdate (cycle_date,jdate,year,jday,hour) !---- Forecast starts from time defined by cycle_date idiff=int(n*dt_sec-dt_sec/2.+0.5) call add_to_jdate(year,jday,hour,idiff,year_fg, & jday_fg,hour_fg) do kk=kk_start(iobs),kk_end(iobs) !----- idiff=(sec), delobs=(min) idiff=kk*delobs*60 !----- ASSUMPTION: observations are defined at equal time intervals (in min) !----- starting from time defined by cycle_date call add_to_jdate(year,jday,hour,idiff,year_obs, & jday_obs,hour_obs) call diff_date1_date2(year_obs,jday_obs,hour_obs, & year_fg,jday_fg,hour_fg,idiff_obs) if(idiff_obs.ge.0.and.idiff_obs.lt.int(dt_sec)) then !----- define output jdate_fg call jdate_make_big (year_obs,jday_obs,hour_obs,jdate_fg) call jdate_to_rdate (jdate_fg,rdate_out) write(file_jdate,1003) OBSTYPE,rdate_out 1003 format('jdate_',a6,'_',a14) CLOSE(i_jdate) OPEN(UNIT=i_jdate,FILE=file_jdate,FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) i_jdate,' OPEN UNIT ERROR IER=',IER read(i_jdate,105,end=107) jdate_obs CLOSE(i_jdate,status='keep') 105 format(a13) if(jdate_fg.eq.jdate_obs) then !------- START model specific part --------- allocate(ensda_work_gbl(1:im,1:jm)) allocate(h_gbl(1:im,1:jm)) allocate(psi_gbl(1:im,1:jm)) allocate(chi_gbl(1:im,1:jm)) call gather(swm_name,im,jm,iim,jjm,1,nsdm,h,ensda_work_gbl) h_gbl(:,:)=ensda_work_gbl(:,:) call gather(swm_name,im,jm,iim,jjm,1,nsdm,psi,ensda_work_gbl) psi_gbl(:,:)=ensda_work_gbl(:,:) call gather(swm_name,im,jm,iim,jjm,1,nsdm,chi,ensda_work_gbl) chi_gbl(:,:)=ensda_work_gbl(:,:) deallocate(ensda_work_gbl) !------- END model specific part --------- if(my_task.eq.0) then write(file_fguess,1001) OBSTYPE,rdate_out 1001 format('fguess_',a6,'.',a14) write(*,*) " file_fguess =",file_fguess CLOSE(i_fguess) OPEN(UNIT=i_fguess,FILE=file_fguess,FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) i_fguess,' OPEN UNIT ERROR IER=',IER do i=1,fgvar_max if(fgvar_name(i).eq.'h ') then write(i_fguess) fgvar_name(i) write(i_fguess) h_gbl end if if(fgvar_name(i).eq.'psi ') then write(i_fguess) fgvar_name(i) write(i_fguess) psi_gbl end if if(fgvar_name(i).eq.'chi ') then write(i_fguess) fgvar_name(i) write(i_fguess) chi_gbl end if end do CLOSE(i_fguess,status='KEEP') end if !! (end my_task=0) deallocate(h_gbl) ; deallocate(psi_gbl) ; deallocate(chi_gbl) else write(*,*) "PROBLEM in write_fg, jdate_fg=",jdate_fg, & "jdate_obs=",jdate_obs endif 107 continue kk_start(iobs)=kk+1 output=.TRUE. exit endif !! if(idiff_obs.ge.0.and.idiff_obs.lt.dt_sec) then if(idiff_obs.gt.0) exit end do !!! kk=kk_start,kobs endif !!! if(OBSFLAG.eq.'YES') then end do !!! iobs=1,iobsmax if(my_task.eq.0) then if(output) then CLOSE(i_parm) OPEN(UNIT=i_parm,FILE='kk_index.parm',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*) i_parm,' OPEN UNIT ERROR IER=',IER write(i_parm,*) iobsmax do i=1,iobsmax write(i_parm,*) kk_start(i),kk_end(i) end do CLOSE(i_parm,status='keep') end if end if !! (end my_task=0) end subroutine write_fg !====================================================================== ! END OF WRITE_FG !====================================================================== !====================================================================== ! BEGINNING OF BIAS_TIME !====================================================================== subroutine bias_time & (n,my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi) ! ********************************************************************** ! ! PROGRAM: bias_time update model bias in current time step ! PRGMMR: M. ZUPANSKI ORG: CIRA/CSU DATE: 2004-05-07 ! ! ! REVISION HISTORY: ! ! 05/07/2004 ..... M. ZUPANSKI ! !----------------------------------------------------------------------- use ensda_variables implicit NONE integer :: my_task integer :: im,jm,iim,jjm,nsdm real(kind=selected_real_kind(13)),dimension(iim,jjm,nsdm) :: h,psi,chi character(*) :: swm_name integer :: i,j,k integer :: n !=========================================================== if(open_bias) then if((n.ne.time).and.(mod(n,nerr).eq.0.or.n.eq.1)) then nbias=nbias+1 call read_model_bias (my_task,swm_name,im,jm,iim,jjm,nsdm) endif !-- Markov process variable (systematic error) if(associated(h_phi_ptr)) then do i=1,iim do j=1,jjm do k=1,nsdm h_phi_ptr(i,j,k)=aa*h_phi_ptr(i,j,k)+(1.-aa)*h_bias_ptr(i,j,k) h(i,j,k)=h(i,j,k)+h_phi_ptr(i,j,k) end do end do end do end if if(associated(psi_phi_ptr)) then do i=1,iim do j=1,jjm do k=1,nsdm psi_phi_ptr(i,j,k)=aa*psi_phi_ptr(i,j,k)+(1.-aa)*psi_bias_ptr(i,j,k) psi(i,j,k)=psi(i,j,k)+psi_phi_ptr(i,j,k) end do end do end do end if if(associated(chi_phi_ptr)) then do i=1,iim do j=1,jjm do k=1,nsdm chi_phi_ptr(i,j,k)=aa*chi_phi_ptr(i,j,k)+(1.-aa)*chi_bias_ptr(i,j,k) chi(i,j,k)=chi(i,j,k)+chi_phi_ptr(i,j,k) end do end do end do end if end if !!!! if(open_bias) then end subroutine bias_time !====================================================================== ! END OF BIAS_TIME !====================================================================== !====================================================================== ! BEGINNING OF READ_MODEL_BIAS !====================================================================== subroutine read_model_bias & (my_task,swm_name,im,jm,iim,jjm,nsdm) ! ********************************************************************** ! ! ROUTINE: read_model_bias: read file with bias ! used in ensda ! PRGMMR: D. ZUPANSKI ORG: CIRA/CSU DATE: 2003-09-17 ! ! ! REVISION HISTORY: ! ! 09/17/2003 ..... D. ZUPANSKI ! !----------------------------------------------------------------------- use gather_scatter use ensda_variables implicit NONE integer :: my_task integer :: im,jm,iim,jjm,nsdm character(*) :: swm_name integer,parameter :: i_bias=410 ! starting bias file integer ubias,unitbfile,IER,i character(len=20) :: file_bias integer :: cvar_max1 character (len=9),dimension(:),allocatable :: cvar_name1 integer,dimension(:),allocatable :: & cvar_num_bias ! number of biases per cntrl vrbl character (len=9),dimension(:),allocatable :: & cvar_name ! cntrl vrbl name logical,dimension(:),allocatable :: & cvar_ic, &! init cond logical mask (per vrbl) cvar_param, &! empir param logical mask (per vrbl) cvar_bias ! model bias logical mask (per vrbl) !======================================================================================== if(open_bias) then !-- read control variable list --------------- CLOSE(103) OPEN(UNIT=103,FILE='cntrl_vrbl',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*)' 103 OPEN UNIT ERROR IER=',IER READ(103,*) cvar_max1 if(cvar_max.ne.cvar_max1) then write(*,*) "PROBLEM in READ_BIAS: cvar_max=",cvar_max," cvar_max1=",cvar_max1 stop else allocate(cvar_name(1:cvar_max)) allocate(cvar_ic(1:cvar_max)) allocate(cvar_param(1:cvar_max)) allocate(cvar_bias(1:cvar_max)) allocate(cvar_num_bias(1:cvar_max)) end if do i=1,cvar_max READ(103,*) cvar_name(i),cvar_ic(i),cvar_param(i),cvar_bias(i),cvar_num_bias(i) end do CLOSE(103,STATUS='KEEP') !------------------------------------------------- ubias=i_bias if(nbias.le.num_bias) then write(file_bias,1003) nbias 1003 format('model_bias_',i2.2) if(my_task.eq.0) then write(*,*) " file_bias =",file_bias end if unitbfile=ubias+nbias CLOSE(unitbfile) OPEN(UNIT=unitbfile,FILE=file_bias,FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) unitbfile,' OPEN UNIT ERROR IER=',IER allocate(ensda_work_gbl(1:im,1:jm)) allocate(ensda_work_loc(1:iim,1:jjm,1:nsdm)) allocate(cvar_name1(1:cvar_max)) do i=1,cvar_max if(cvar_bias(i)) then read(unitbfile) cvar_name1(i) if(cvar_name(i).ne.cvar_name1(i)) then write(*,*) "PROBLEM in READ_BIAS: cvar_name,cvar_name1=",cvar_name(i),cvar_name1(i) stop end if if(cvar_name(i).eq.'h ') then allocate(h_bias_gbl(1:im,1:jm)) read(unitbfile) h_bias_gbl ensda_work_gbl(:,:)=h_bias_gbl(:,:) call scatter(swm_name,im,jm,iim,jjm,1,nsdm,ensda_work_gbl,ensda_work_loc) h_bias(:,:,:)=ensda_work_loc(:,:,:) deallocate(h_bias_gbl) end if if(cvar_name(i).eq.'psi ') then allocate(psi_bias_gbl(1:im,1:jm)) read(unitbfile) psi_bias_gbl ensda_work_gbl(:,:)=psi_bias_gbl(:,:) call scatter(swm_name,im,jm,iim,jjm,1,nsdm,ensda_work_gbl,ensda_work_loc) psi_bias(:,:,:)=ensda_work_loc(:,:,:) deallocate(psi_bias_gbl) end if if(cvar_name(i).eq.'chi ') then allocate(chi_bias_gbl(1:im,1:jm)) read(unitbfile) chi_bias_gbl ensda_work_gbl(:,:)=chi_bias_gbl(:,:) call scatter(swm_name,im,jm,iim,jjm,1,nsdm,ensda_work_gbl,ensda_work_loc) chi_bias(:,:,:)=ensda_work_loc(:,:,:) deallocate(chi_bias_gbl) end if end if end do !!! do i=1,cvar_max CLOSE(unitbfile,status='KEEP') deallocate(cvar_name1) deallocate(ensda_work_gbl) deallocate(ensda_work_loc) else if(my_task.eq.0) then write(*,*)"read_model_bias PROBLEM nbias=",nbias," num_bias=",num_bias end if endif !!! if(nbias.le.num_bias) deallocate(cvar_name) deallocate(cvar_ic) deallocate(cvar_param) deallocate(cvar_bias) deallocate(cvar_num_bias) end if !!! if(open_bias) end subroutine read_model_bias !====================================================================== ! END OF READ_MODEL_BIAS !====================================================================== !====================================================================== ! BEGINNING OF READ_MODEL_PARAM !====================================================================== subroutine read_model_param & (my_task,visc) ! ********************************************************************** ! ! ROUTINE: read_model_param: read file with param ! used in ensda ! PRGMMR: D. ZUPANSKI ORG: CIRA/CSU DATE: 2003-09-17 ! ! ! REVISION HISTORY: ! ! 09/17/2003 ..... D. ZUPANSKI ! !----------------------------------------------------------------------- use ensda_variables implicit NONE integer :: my_task real :: visc integer,parameter :: i_param=403 ! empirical parameter file integer :: IER,i integer :: cvar_max1 character (len=9),dimension(:),allocatable :: cvar_name1 integer,dimension(:),allocatable :: & cvar_num_bias ! number of biases per cntrl vrbl character (len=9),dimension(:),allocatable :: & cvar_name ! cntrl vrbl name logical,dimension(:),allocatable :: & cvar_ic, &! init cond logical mask (per vrbl) cvar_param, &! empir param logical mask (per vrbl) cvar_bias ! model bias logical mask (per vrbl) !======================================================================================= if(open_param) then !-- read control variable list --------------- CLOSE(103) OPEN(UNIT=103,FILE='cntrl_vrbl',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*)' 103 OPEN UNIT ERROR IER=',IER READ(103,*) cvar_max1 if(cvar_max.ne.cvar_max1) then write(*,*) "PROBLEM in READ_PARAM: cvar_max=",cvar_max," cvar_max1=",cvar_max1 stop else allocate(cvar_name(1:cvar_max)) allocate(cvar_ic(1:cvar_max)) allocate(cvar_param(1:cvar_max)) allocate(cvar_bias(1:cvar_max)) allocate(cvar_num_bias(1:cvar_max)) end if do i=1,cvar_max READ(103,*) cvar_name(i),cvar_ic(i),cvar_param(i),cvar_bias(i),cvar_num_bias(i) end do CLOSE(103,STATUS='KEEP') !------------------------------------------------- CLOSE(i_param) OPEN(UNIT=i_param,FILE='model_param',FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) i_param,' OPEN UNIT ERROR IER=',IER allocate(cvar_name1(1:cvar_max)) do i=1,cvar_max if(cvar_param(i)) then read(i_param) cvar_name1(i) if(cvar_name(i).ne.cvar_name1(i)) then write(*,*) "PROBLEM in READ_PARAM: cvar_name,cvar_name1=",cvar_name(i),cvar_name1(i) stop end if if(cvar_name(i).eq.'visc ') then read(i_param) visc if(my_task.eq.0) then if(visc.lt.0.) then write(*,*) "NEGATIVE VISCOUS PARAMETER visc=",visc write(*,*) "CONTINUE WITH VISCOUS PARAMETER visc=0.0" visc=0.0 else write(*,*) "CURRENT VISCOUS PARAMETER visc=",visc end if end if end if end if enddo deallocate(cvar_name1) CLOSE(i_param,status='KEEP') deallocate(cvar_name) deallocate(cvar_ic) deallocate(cvar_param) deallocate(cvar_bias) deallocate(cvar_num_bias) end if !!! if(open_param) end subroutine read_model_param !====================================================================== ! END OF READ_MODEL_PARAM !====================================================================== !====================================================================== ! BEGINNING OF READ_MODEL_IC !====================================================================== subroutine read_model_ic & (my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi) ! ********************************************************************** ! ! ROUTINE: read_model_ic: read file with ic ! used in ensda ! PRGMMR: D. ZUPANSKI ORG: CIRA/CSU DATE: 2003-09-17 ! ! ! REVISION HISTORY: ! ! 09/17/2003 ..... D. ZUPANSKI ! !----------------------------------------------------------------------- use gather_scatter use ensda_variables implicit NONE integer :: my_task integer :: im,jm,iim,jjm,nsdm real(kind=selected_real_kind(13)),dimension(iim,jjm,nsdm) :: h,psi,chi character(*) :: swm_name integer,parameter :: i_ic=402 ! initial conditions file integer IER,i integer :: cvar_max1 character (len=9),dimension(:),allocatable :: cvar_name1 integer,dimension(:),allocatable :: & cvar_num_bias ! number of biases per cntrl vrbl character (len=9),dimension(:),allocatable :: & cvar_name ! cntrl vrbl name logical,dimension(:),allocatable :: & cvar_ic, &! init cond logical mask (per vrbl) cvar_param, &! empir param logical mask (per vrbl) cvar_bias ! model bias logical mask (per vrbl) !============================================================================= if(open_ic) then !-- read control variable list --------------- CLOSE(103) OPEN(UNIT=103,FILE='cntrl_vrbl',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*)' 103 OPEN UNIT ERROR IER=',IER READ(103,*) cvar_max1 if(cvar_max.ne.cvar_max1) then write(*,*) "PROBLEM in READ_IC: cvar_max=",cvar_max," cvar_max1=",cvar_max1 stop else allocate(cvar_name(1:cvar_max)) allocate(cvar_ic(1:cvar_max)) allocate(cvar_param(1:cvar_max)) allocate(cvar_bias(1:cvar_max)) allocate(cvar_num_bias(1:cvar_max)) end if do i=1,cvar_max READ(103,*) cvar_name(i),cvar_ic(i),cvar_param(i),cvar_bias(i),cvar_num_bias(i) end do CLOSE(103,STATUS='KEEP') !------------------------------------------------- CLOSE(i_ic) OPEN(UNIT=i_ic,FILE='model_ic',FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) i_ic,' OPEN UNIT ERROR IER=',IER allocate(ensda_work_gbl(1:im,1:jm)) allocate(cvar_name1(1:cvar_max)) do i=1,cvar_max if(cvar_ic(i)) then read(i_ic) cvar_name1(i) if(cvar_name(i).ne.cvar_name1(i)) then write(*,*) "PROBLEM in READ_IC: cvar_name,cvar_name1=",cvar_name(i),cvar_name1(i) stop end if if(cvar_name(i).eq.'h ') then allocate(h_gbl(1:im,1:jm)) read(i_ic) h_gbl ensda_work_gbl(:,:)=h_gbl(:,:) call scatter(swm_name,im,jm,iim,jjm,1,nsdm,ensda_work_gbl,h) deallocate(h_gbl) end if if(cvar_name(i).eq.'psi ') then allocate(psi_gbl(1:im,1:jm)) read(i_ic) psi_gbl ensda_work_gbl(:,:)=psi_gbl(:,:) call scatter(swm_name,im,jm,iim,jjm,1,nsdm,ensda_work_gbl,psi) deallocate(psi_gbl) end if if(cvar_name(i).eq.'chi ') then allocate(chi_gbl(1:im,1:jm)) read(i_ic) chi_gbl ensda_work_gbl(:,:)=chi_gbl(:,:) call scatter(swm_name,im,jm,iim,jjm,1,nsdm,ensda_work_gbl,chi) deallocate(chi_gbl) end if end if end do CLOSE(i_ic,status='KEEP') deallocate(ensda_work_gbl) deallocate(cvar_name1) deallocate(cvar_name) deallocate(cvar_ic) deallocate(cvar_param) deallocate(cvar_bias) deallocate(cvar_num_bias) end if !!! if(open_ic) end subroutine read_model_ic !====================================================================== ! END OF READ_MODEL_IC !====================================================================== !====================================================================== ! BEGINNING OF WRITE_MODEL_BIAS !====================================================================== subroutine write_model_bias & (my_task,swm_name,im,jm,iim,jjm,nsdm) ! ********************************************************************** ! ! ROUTINE: write_model_bias: write file with bias ! used in ensda ! PRGMMR: D. ZUPANSKI ORG: CIRA/CSU DATE: 2003-09-17 ! ! ! REVISION HISTORY: ! ! 09/17/2003 ..... D. ZUPANSKI ! !----------------------------------------------------------------------- use gather_scatter use ensda_variables implicit NONE integer :: my_task integer :: im,jm,iim,jjm,nsdm character(*) :: swm_name integer,parameter :: i_bias=410 ! starting bias file integer ubias,unitbfile,IER,i character(len=20) :: file_bias integer :: cvar_max1 integer,dimension(:),allocatable :: & cvar_num_bias ! number of biases per cntrl vrbl character (len=9),dimension(:),allocatable :: & cvar_name ! cntrl vrbl name logical,dimension(:),allocatable :: & cvar_ic, &! init cond logical mask (per vrbl) cvar_param, &! empir param logical mask (per vrbl) cvar_bias ! model bias logical mask (per vrbl) !================================================================================ if(open_bias) then !-- read control variable list --------------- CLOSE(103) OPEN(UNIT=103,FILE='cntrl_vrbl',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*)' 103 OPEN UNIT ERROR IER=',IER READ(103,*) cvar_max1 if(cvar_max.ne.cvar_max1) then write(*,*) "PROBLEM in WRITE_BIAS: cvar_max=",cvar_max," cvar_max1=",cvar_max1 stop else allocate(cvar_name(1:cvar_max)) allocate(cvar_ic(1:cvar_max)) allocate(cvar_param(1:cvar_max)) allocate(cvar_bias(1:cvar_max)) allocate(cvar_num_bias(1:cvar_max)) end if do i=1,cvar_max READ(103,*) cvar_name(i),cvar_ic(i),cvar_param(i),cvar_bias(i),cvar_num_bias(i) end do CLOSE(103,STATUS='KEEP') !------------------------------------------------- ubias=i_bias if(nbias.le.num_bias) then write(file_bias,1003) nbias 1003 format('model_bias_',i2.2) if(my_task.eq.0) then write(*,*) " file_bias =",file_bias end if unitbfile=ubias+nbias if(my_task.eq.0) then CLOSE(unitbfile) OPEN(UNIT=unitbfile,FILE=file_bias,FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) unitbfile,' OPEN UNIT ERROR IER=',IER end if allocate(ensda_work_gbl(1:im,1:jm)) allocate(ensda_work_loc(1:iim,1:jjm,1:nsdm)) do i=1,cvar_max if(cvar_bias(i)) then if(cvar_name(i).eq.'h ') then allocate(h_bias_gbl(1:im,1:jm)) ensda_work_loc(:,:,:)=h_bias(:,:,:) call gather(swm_name,im,jm,iim,jjm,1,nsdm,ensda_work_loc,ensda_work_gbl) h_bias_gbl(:,:)=ensda_work_gbl(:,:) if(my_task.eq.0) then write(unitbfile) cvar_name(i) write(unitbfile) h_bias_gbl end if deallocate(h_bias_gbl) endif if(cvar_name(i).eq.'psi ') then allocate(psi_bias_gbl(1:im,1:jm)) ensda_work_loc(:,:,:)=psi_bias(:,:,:) call gather(swm_name,im,jm,iim,jjm,1,nsdm,ensda_work_loc,ensda_work_gbl) psi_bias_gbl(:,:)=ensda_work_gbl(:,:) if(my_task.eq.0) then write(unitbfile) cvar_name(i) write(unitbfile) psi_bias_gbl end if deallocate(psi_bias_gbl) endif if(cvar_name(i).eq.'chi ') then allocate(chi_bias_gbl(1:im,1:jm)) ensda_work_loc(:,:,:)=chi_bias(:,:,:) call gather(swm_name,im,jm,iim,jjm,1,nsdm,ensda_work_loc,ensda_work_gbl) chi_bias_gbl(:,:)=ensda_work_gbl(:,:) if(my_task.eq.0) then write(unitbfile) cvar_name(i) write(unitbfile) chi_bias_gbl end if deallocate(chi_bias_gbl) endif endif end do !!! do i=1,cvar_max if(my_task.eq.0) then CLOSE(unitbfile,status='KEEP') end if deallocate(ensda_work_gbl) deallocate(ensda_work_loc) else if(my_task.eq.0) then write(*,*)"write_model_bias PROBLEM nbias=",nbias," num_bias=",num_bias end if endif !!! if(nbias.le.num_bias) deallocate(cvar_name) deallocate(cvar_ic) deallocate(cvar_param) deallocate(cvar_bias) deallocate(cvar_num_bias) end if !!! if(open_bias) end subroutine write_model_bias !====================================================================== ! END OF WRITE_MODEL_BIAS !====================================================================== !====================================================================== ! BEGINNING OF WRITE_MODEL_PARAM !====================================================================== subroutine write_model_param & (my_task,visc) ! ********************************************************************** ! ! ROUTINE: write_model_param: write file with param ! used in ensda ! PRGMMR: D. ZUPANSKI ORG: CIRA/CSU DATE: 2003-09-17 ! ! ! REVISION HISTORY: ! ! 09/17/2003 ..... D. ZUPANSKI ! !----------------------------------------------------------------------- use ensda_variables implicit NONE integer :: my_task real :: visc integer,parameter :: i_param=403 ! empirical parameter file integer IER,i integer :: cvar_max1 integer,dimension(:),allocatable :: & cvar_num_bias ! number of biases per cntrl vrbl character (len=9),dimension(:),allocatable :: & cvar_name ! cntrl vrbl name logical,dimension(:),allocatable :: & cvar_ic, &! init cond logical mask (per vrbl) cvar_param, &! empir param logical mask (per vrbl) cvar_bias ! model bias logical mask (per vrbl) !=========================================================================== if(open_param) then !-- read control variable list --------------- CLOSE(103) OPEN(UNIT=103,FILE='cntrl_vrbl',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*)' 103 OPEN UNIT ERROR IER=',IER READ(103,*) cvar_max1 if(cvar_max.ne.cvar_max1) then write(*,*) "PROBLEM in WRITE_PARAM: cvar_max=",cvar_max," cvar_max1=",cvar_max1 stop else allocate(cvar_name(1:cvar_max)) allocate(cvar_ic(1:cvar_max)) allocate(cvar_param(1:cvar_max)) allocate(cvar_bias(1:cvar_max)) allocate(cvar_num_bias(1:cvar_max)) end if do i=1,cvar_max READ(103,*) cvar_name(i),cvar_ic(i),cvar_param(i),cvar_bias(i),cvar_num_bias(i) end do CLOSE(103,STATUS='KEEP') !------------------------------------------------- if(my_task.eq.0) then CLOSE(i_param) OPEN(UNIT=i_param,FILE='model_param',FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) i_param,' OPEN UNIT ERROR IER=',IER end if do i=1,cvar_max if(cvar_param(i)) then if(cvar_name(i).eq.'visc ') then if(my_task.eq.0) then write(i_param) cvar_name(i) write(i_param) visc end if end if end if enddo if(my_task.eq.0) then CLOSE(i_param,status='KEEP') end if deallocate(cvar_name) deallocate(cvar_ic) deallocate(cvar_param) deallocate(cvar_bias) deallocate(cvar_num_bias) end if !!! if(open_param) end subroutine write_model_param !====================================================================== ! END OF WRITE_MODEL_PARAM !====================================================================== !====================================================================== ! BEGINNING OF WRITE_MODEL_IC !====================================================================== subroutine write_model_ic & (my_task,swm_name,im,jm,iim,jjm,nsdm,h,psi,chi) ! ********************************************************************** ! ! ROUTINE: write_model_ic: write file with ic ! used in ensda ! PRGMMR: D. ZUPANSKI ORG: CIRA/CSU DATE: 2003-09-17 ! ! ! REVISION HISTORY: ! ! 09/17/2003 ..... D. ZUPANSKI ! !----------------------------------------------------------------------- use gather_scatter use ensda_variables implicit NONE integer :: my_task integer :: im,jm,iim,jjm,nsdm real(kind=selected_real_kind(13)),dimension(iim,jjm,nsdm) :: h,psi,chi character(*) :: swm_name integer,parameter :: i_ic=402 ! initial conditions file integer IER,i integer i1,i2 integer :: cvar_max1 integer,dimension(:),allocatable :: & cvar_num_bias ! number of biases per cntrl vrbl character(len=9),dimension(:),allocatable :: & cvar_name ! cntrl vrbl name logical,dimension(:),allocatable :: & cvar_ic, &! init cond logical mask (per vrbl) cvar_param, &! empir param logical mask (per vrbl) cvar_bias ! model bias logical mask (per vrbl) !=============================================================================== if(open_ic) then !-- read control variable list --------------- CLOSE(103) OPEN(UNIT=103,FILE='cntrl_vrbl',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*)' 103 OPEN UNIT ERROR IER=',IER READ(103,*) cvar_max1 if(cvar_max.ne.cvar_max1) then write(*,*) "PROBLEM in WRITE_IC: cvar_max=",cvar_max," cvar_max1=",cvar_max1 stop else allocate(cvar_name(1:cvar_max)) allocate(cvar_ic(1:cvar_max)) allocate(cvar_param(1:cvar_max)) allocate(cvar_bias(1:cvar_max)) allocate(cvar_num_bias(1:cvar_max)) end if do i=1,cvar_max READ(103,*) cvar_name(i),cvar_ic(i),cvar_param(i),cvar_bias(i),cvar_num_bias(i) end do CLOSE(103,STATUS='KEEP') !------------------------------------------------- if(my_task.eq.0) then CLOSE(i_ic) OPEN(UNIT=i_ic,FILE='model_ic',FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) i_ic,' OPEN UNIT ERROR IER=',IER end if allocate(ensda_work_gbl(1:im,1:jm)) do i=1,cvar_max if(cvar_ic(i)) then if(cvar_name(i).eq.'h ') then allocate(h_gbl(1:im,1:jm)) call gather(swm_name,im,jm,iim,jjm,1,nsdm,h,ensda_work_gbl) h_gbl(:,:)=ensda_work_gbl(:,:) if(my_task.eq.0) then write(i_ic) cvar_name(i) write(i_ic) ((h_gbl(i1,i2),i1=1,im),i2=1,jm) end if deallocate(h_gbl) end if if(cvar_name(i).eq.'psi ') then allocate(psi_gbl(1:im,1:jm)) call gather(swm_name,im,jm,iim,jjm,1,nsdm,psi,ensda_work_gbl) psi_gbl(:,:)=ensda_work_gbl(:,:) if(my_task.eq.0) then write(i_ic) cvar_name(i) write(i_ic) ((psi_gbl(i1,i2),i1=1,im),i2=1,jm) end if deallocate(psi_gbl) end if if(cvar_name(i).eq.'chi ') then allocate(chi_gbl(1:im,1:jm)) call gather(swm_name,im,jm,iim,jjm,1,nsdm,chi,ensda_work_gbl) chi_gbl(:,:)=ensda_work_gbl(:,:) if(my_task.eq.0) then write(i_ic) cvar_name(i) write(i_ic) ((chi_gbl(i1,i2),i1=1,im),i2=1,jm) end if deallocate(chi_gbl) end if end if end do if(my_task.eq.0) then CLOSE(i_ic,status='KEEP') end if deallocate(ensda_work_gbl) deallocate(cvar_name) deallocate(cvar_ic) deallocate(cvar_param) deallocate(cvar_bias) deallocate(cvar_num_bias) end if !!! if(open_ic) end subroutine write_model_ic !====================================================================== ! END OF WRITE_MODEL_IC !======================================================================