PROGRAM Rcm USE Rcm_mod_subs IMPLICIT NONE ! CHARACTER(LEN=8) :: real_date CHARACTER (LEN=8) :: time_char CHARACTER(LEN=10) ::real_time CHARACTER(LEN=80) :: ST='', PS='', HD='' LOGICAL :: FD ! ! INTEGER (iprec) :: k, kc, n REAL (rprec) :: dt ! ! IF (.NOT.Check_logical_units ( ) ) STOP 'LUNs NOT AVAILABLE' ! ! ! Grid limit parameters i1 = 2 ! will reset anyway i2 = isize - 1 j1 = jwrap j2 = jsize - 1 iint = 1 jint = 1 ! ! ! !**** Read input data sets: ****** ! CALL Read_grid () ! CALL Read_plasma () ! CALL Rcm_read_params () ! ! Read in other inputs, both constant and run-specific: ! CALL Read_vdrop () CALL Read_kp () CALL Read_bfield () print*,'bfield' CALL Read_eta_on_bndy () CALL Read_qtcond () CALL Read_winds (iwind) CALL Read_dktime (L_dktime) CALL Read_trf () ! ! ! ! ! Open file for formatted output and do initial print out : ! CALL Initial_printout () ! ! ! i1 = imin + 1 ! ! !--> Read initial locations of inner edges ! IF (itimei == 0) THEN CALL Read_array ('rcmbi_inp', irdr, label, ARRAY_1D = bi, ASCI = asci_flag) CALL Read_array ('rcmbj_inp', irdr, label, ARRAY_1D = bj, ASCI = asci_flag) CALL Read_array ('rcmetab_inp', irdr, label, ARRAY_1D = etab, ASCI = asci_flag) CALL Read_array ('rcmitrack_inp', irdr, label, ARRAY_1D = itrack, ASCI = asci_flag) CALL Read_array ('rcmmpoint_inp', irdr, label, ARRAY_1D = mpoint, ASCI = asci_flag) CALL Read_array ('rcmnpoint_inp', irdr, label, ARRAY_1D = npoint, ASCI = asci_flag) CALL Read_array ('rcmeeta_inp', irdr, label, ARRAY_3D = eeta, ASCI = asci_flag) ELSE CALL Read_array ('rcmbi', irdr, label, ARRAY_1D = bi) CALL Read_array ('rcmbj', irdr, label, ARRAY_1D = bj) CALL Read_array ('rcmetab', irdr, label, ARRAY_1D = etab) CALL Read_array ('rcmitrack', irdr, label, ARRAY_1D = itrack) CALL Read_array ('rcmmpoint', irdr, label, ARRAY_1D = mpoint) CALL Read_array ('rcmnpoint', irdr, label, ARRAY_1D = npoint) CALL Read_array ('rcmeeta', irdr, label, ARRAY_3D = eeta) END IF ! ! itrack (nptmax) = nptmax ! 3/96 frt IF (MAXVAL (ABS(itrack)) > nptmax) THEN WRITE (*,*) 'RCM: itrack(nptmax) is being reset' itrack (nptmax) = MAXVAL (ABS(itrack)) + 1 END IF ! ! ! ! IF hot restart, read V and check the time label: ! IF (itimei /= 0) THEN WRITE (*,'(A)', ADVANCE='NO') & 'HOT restart, reading V from file to check time label...' CALL READ_array ('rcmv', irdr, label, ARRAY_2D = v) IF (label%intg(6) /= itimei ) STOP 'T in file/=ITIMEI V' WRITE (*,*) 'OK' !$ WRITE (*,*) 'RESETTING V TO ZERO'; v = 0.0 ! END IF ! ! ! ! !******************* main time loop ************************* ! IF (idt1/idt*idt /= idt1) STOP 'RCM: idt1--idt' itout1 = itimei ! next time to write to disk; set this to ! ITIMEI to write out initial configuration itout2 = itimei ! next time to do formatted output itcln = itimei ! next time to call ADD & ZAP ! dt = REAL (idt) ! ! !: print*,'itimei=',itimei,itimef,idt DO i_time = itimei, itimef, idt ! print*,'i_time=',i_time ! DO j_time = 1, idt_reduce_factor dt = REAL (idt) / REAL (idt_reduce_factor) ! time = itime + idt/idt_reduce_factor*(j_time-1) ! CALL Comput (i_time, dt) ! !$ CALL Disk_write_arrays () IF (i_time == itout1 .OR. i_time == itimef) THEN CALL Disk_write_arrays () IF (irdw > 1 .AND. L_spiro) CALL Spiro_effect (i_time, idt, irdw) itout1 = itout1 + idt1 irdw = irdw + 1 END IF ! ! ! sazykin 3/25/05: don't touch edges if they are not used! IF (i_time == itcln) THEN IF (L_move_plasma_edges) CALL Clean_up_edges () itcln = itcln + idt3 END IF ! ! IF(i_time==itout2 .OR. i_time==itimef) THEN CALL Formatted_output () END IF ! ! CALL Move_plasma ( dt ) ! END DO END DO ! ! CLOSE (UNIT = LUN_2) CLOSE (UNIT = LUN_3) STOP ! ! CONTAINS ! ! SUBROUTINE Initial_printout () INTEGER (iprec) :: j ! CALL Date_and_time (real_date, real_time) IF (itimei == 0) THEN ST = 'REPLACE' PS = 'REWIND' HD = 'BEGINNING NEW RUN' ELSE ST = 'OLD' PS = 'APPEND' HD = 'CONTINUE SAME RUN' END IF OPEN (LUN_2, FILE = 'rcm.printout', STATUS = ST, POSITION = PS) OPEN (LUN_3, FILE = 'rcm.index', STATUS = ST, POSITION = PS) WRITE (LUN_3,'(T2,A)',ADVANCE='NO') TRIM(HD) WRITE (LUN_3,'(A11,A4,A1,A2,A1,A2, A8,A2,A1,A2,A1,A2)') & ' TODAY IS ', real_date(1:4), '/', & real_date(5:6), '/', & real_date (7:8), & ' TIME: ', real_time(1:2), ':', & real_time(3:4), ':', & real_time(5:6) WRITE (LUN_3,902) ! WRITE (LUN_2,*) 'START OF RCM RUN:' WRITE (LUN_2,'(A,I6,A,I6)') 'WILL START AT ITIMEI=',itimei, & ' AND STOP AT ITIMEF=',itimef WRITE (LUN_2,'(T5,A,T35,I5.5)') 'time step =', idt WRITE (LUN_2,'(T5,A,T35,I5.5)') 'disk write time step=', idt1 WRITE (LUN_2,'(T5,A,T35,I5.5)') 'printout time step=', idt2 WRITE (LUN_2,'(T5,A,T35,I5.5)') 'edge cleanup time step=', idt3 WRITE (LUN_2,'(T5,A,T35,I5.5)') 'imin =', imin WRITE (LUN_2,'(T5,A,T35,I5.5)') 'start at itimei=', itimei WRITE (LUN_2,'(T5,A,T35,I5.5)') 'stop at itimef=', itimef WRITE (LUN_2,'(T5,A,T35,I5.5)') 'read at itimei from REC ',irdr WRITE (LUN_2,'(T5,A,T35,I5.5)') 'start writing at REC ', irdw WRITE (LUN_2,'(/T5,A)' ) 'SIZES PARAMETERS:' WRITE (LUN_2,'(T10,A,T20,I6)') 'isize=',isize WRITE (LUN_2,'(T10,A,T20,I6)') 'jsize=',jsize WRITE (LUN_2,'(T10,A,T20,I6)') 'ksize=',ksize WRITE (LUN_2,'(T10,A,T20,I6)') 'kcsize=',kcsize WRITE (LUN_2,'(T10,A,T20,I6)') 'iesize=',iesize WRITE (LUN_2,'(/T5,A)' ) 'GRID PARAMETERS:' WRITE (LUN_2,'(T10,A,T20,G9.2)') 'dlam=',dlam WRITE (LUN_2,'(T10,A,T20,G9.2)') 'dpsi=',dpsi WRITE (LUN_2,'(T10,A,T20,G9.2)') 're=',re WRITE (LUN_2,'(T10,A,T20,G9.2)') 'ri=',ri WRITE (LUN_2,'(T10,A,T20,2G9.2)') 'xmass',xmass WRITE (LUN_2,'(/T5,A)' ) 'PLASMA EDGES PARAMETERS:' DO k = 1, ksize WRITE (LUN_2,'(T10,A,I3,T25,A,I3,T40,A,G9.2,T60,A,ES9.2)')& 'k=', k, 'ikflav=',ikflav(k), & 'alam=', alam(k), 'eta=', eta(k) END DO WRITE (LUN_2,'(/T5,A)' ) 'PLASMA GRID PARAMETERS:' DO kc = 1, kcsize DO j = 1, jsize WRITE (LUN_2,'(T10,A,I3,T20,A,G9.2,T45,A,ES9.2, T65,A,F5.3)') & 'kc=', kc, 'alamc=', alamc(kc), 'etac=', etac(kc,j), 'f=', fudgec(kc) END DO END DO WRITE (LUN_2,'(/T5,A)' ) 'PCP INPUT PARAMETERS:' WRITE (LUN_2,'(T5,A,T20,I5)') 'nvmax =', SIZE(ivtime) DO n = 1, SIZE(ivtime) WRITE (LUN_2,'(T5,A,I7,T20,A,G9.2)') & 'ivtime=',ivtime(n), 'vinput=',vinput(n) END DO WRITE (LUN_2,'(/T5,A)' ) 'BFIELD MARKTIMES:' WRITE (LUN_2,'(T5,A,T20,I6)') 'nbf =', SIZE(ibtime) DO n = 1, SIZE(ibtime) WRITE (LUN_2,'(T5,A,I7)') 'ibtime=',ibtime(n) END DO WRITE (LUN_2,'(/T5,A)' ) 'KP MARKTIMES:' WRITE (LUN_2,'(T5,A,T20,I6)') 'nkpmax =', SIZE(ikptime) DO n = 1, SIZE(ikptime) WRITE (LUN_2,'(T5,A,I7)') 'ikptime=',ikptime(n) END DO WRITE (LUN_2,'(/T5,A)' ) 'Plasma MARKTIMES:' WRITE (LUN_2,'(T5,A,T20,I6)') 'itime_etac =', SIZE(itime_etac) DO n = 1, SIZE(itime_etac) WRITE (LUN_2,'(T5,A,I7)') 'itime_etac=',itime_etac(n) END DO WRITE (LUN_2,'(T2,A,T20,I2)') 'IPOT =', ipot WRITE (LUN_2,'(T2,A,T20,I2)') 'ICOND = ', icond WRITE (LUN_2,'(T2,A,T20,I2)') 'IBND = ', ibnd_type WRITE (LUN_2,'(T2,A,T20,I2)') 'IPCP_TYPE=', ipcp_type ! ! ! 902 FORMAT (T2,'TIME', T12,'ITIME' , T19,'REC#' ,& T26,'VDROP', T33,'KP', T39,'FSTOFF', & T46,'FMEB', T53, 'DST', T62,'FCLPS', T69,'VDROP_PHASE' ) RETURN END SUBROUTINE Initial_printout ! ! ! SUBROUTINE Disk_write_arrays () ! ! Writing rcm arrays to files is done at each time step, but the ! record number is changed only with the time step specified in ! 'rcm.params'. This ensures that if the model crashes, files ! contain the most recent arrays. ! ! !________We call OUTPUT subroutine with this flag. The policy is ! that if we start at time=0, then we delete any old files ! and start from scratch. Otherwise, continue to output to ! existing files if they exist or create them if not. ! IF (i_time == 0) THEN FD = .TRUE. ELSE FD = .FALSE. END IF ! CALL Rcm_write_label () ! WRITE (time_char,'(I2.2,A1,I2.2,A1,I2.2)') & label%intg(3), ':', label%intg(4), ':', label%intg(5) WRITE (*,'(T2,A21,I6.6,A10,TR1,F5.3,A1)') & '-->TIME_STEP, T=', i_time,'('//time_char,& label%real(05),')' ! IF (i_time == itout1 .OR. i_time == itimef) THEN ! close (lun_3) ! open (lun_3, status='old', position='append') WRITE (lun_3,901) time_char, & i_time, irdw, vdrop, kp, fstoff, fmeb, fdst, fclps, vdrop_phase 901 FORMAT (T2,A8, T12,I6, T19,I5, T26,F5.1, T33,F4.2,& T39,F5.2, T46,F5.2, T53,F7.1, T62,F5.1, T69, F6.2) CALL Flush (lun_3) END IF ! IF (idebug == 0) THEN ! ! ! 1. Write magnetic field and bndy location: ! CALL Write_array ('rcmxmin', irdw, label, ARRAY_2D = xmin, SETUP = FD) CALL Write_array ('rcmymin', irdw, label, ARRAY_2D = ymin, SETUP = FD) CALL Write_array ('rcmvm', irdw, label, ARRAY_2D = vm, SETUP = FD) CALL Write_array ('rcmbmin', irdw, label, ARRAY_2D = bmin, SETUP = FD) CALL Write_array ('rcmbndloc', irdw, label, ARRAY_1D = bndloc, SETUP = FD) ! ! ! 2. Write out plasma info: ! CALL Write_array ('rcmetac', irdw, label, ARRAY_2D = etac, SETUP=FD) CALL Write_array ('rcmbi', irdw, label, ARRAY_1D = bi, SETUP=FD ) CALL Write_array ('rcmbj', irdw, label, ARRAY_1D = bj , SETUP=FD) CALL Write_array ('rcmetab', irdw, label, ARRAY_1D = etab , SETUP=FD) CALL Write_array ('rcmitrack', irdw, label, ARRAY_1D = itrack, SETUP=FD ) CALL Write_array ('rcmmpoint', irdw, label, ARRAY_1D = mpoint, SETUP=FD ) CALL Write_array ('rcmnpoint', irdw, label, ARRAY_1D = npoint, SETUP=FD ) CALL Write_array ('rcmeeta', irdw, label, ARRAY_3D = eeta, SETUP=FD ) print*, 'rcmf90eta=',eeta(10,27,2), eeta(20,39,40) ! CALL Write_array ('rcmdbidt', irdw, label, ARRAY_1D = dbidt, SETUP = FD ) ! CALL Write_array ('rcmdbjdt', irdw, label, ARRAY_1D = dbjdt, SETUP = FD ) ! ! ! ! 3. Write ionospheric quantities: ! CALL Write_array ('rcmpedlam', irdw, label, ARRAY_2D = pedlam, SETUP=FD ) CALL Write_array ('rcmpedpsi', irdw, label, ARRAY_2D = pedpsi, SETUP=FD ) CALL Write_array ('rcmhall', irdw, label, ARRAY_2D = hall , SETUP=FD ) CALL Write_array ('rcmssequat',irdw, label, ARRAY_1D = ss , SETUP=FD ) CALL Write_array ('rcmeavg', irdw, label, ARRAY_3D = eavg , SETUP=FD ) CALL Write_array ('rcmeflux', irdw, label, ARRAY_3D = eflux , SETUP=FD ) CALL Write_array ('rcmbirk', irdw, label, ARRAY_2D = birk , SETUP=FD ) ! ! ! 4. Write out magnetospheric quantities: ! CALL Write_array ('rcmv', irdw, label, ARRAY_2D = v, SETUP=FD ) CALL Write_array ('rcmvpar', irdw, label, ARRAY_2D = vpar, SETUP=FD ) ! END IF RETURN END SUBROUTINE Disk_write_arrays ! ! ! SUBROUTINE Clean_up_edges () IMPLICIT NONE INTEGER (iprec) :: k ! ! Check if pts need to be deleted or added: ! DO k = 1, ksize ! CALL Farmers_wife (k) ! CALL Zap_edge (k) ! CALL Expand_edge (k) ! !!! sts CALL set_itrack (jdim, kdim, nptmax, bi, bj, itrack, mpoint, & !!! npoint,ain) ! END DO RETURN END SUBROUTINE Clean_up_edges ! ! ! SUBROUTINE Formatted_output () WRITE (LUN_2, *) irdw, i_time itout2 = itout2 + idt2 RETURN END SUBROUTINE Formatted_output ! ! END PROGRAM Rcm