#define ITYPE_RST 0 #define MACHINE_WORD 4 #define c_str(s) ('s\0') subroutine init_rstrt (nxp, nyp, nz, npt, zin, iseg) c--------------------------------------------------------- include 'comm_new.h' include 'comm_data.h' character*100 str common /run/ nstart,nlaststart,nskip,nsteps,nergy,nskipo,nlast common /new_save/ iword, nruns, nscpu, nswll, inf(100), rnf(100) dimension zin(1), zz(50), iseg(1) equivalence (inf(20),kd_sta), * (inf(21),kd_xy), (inf(22),kd_z), (inf(23),kd_time), * (inf(24),kd_seg), (inf(25),kd_uv), (inf(26),kd_temp), * (inf(27),kd_salt),(inf(28),kd_tr), (inf(29),kd_phi), * (inf(30),kd_conv), (inf(31),kd_means), * (inf(100),kd_end) call opda (1, 1, fbo(1:n_out)) do i = 1, 100 inf(i) = 0 rnf(i) = 0. enddo iword = MACHINE_WORD kd = iword * 200 str = 'Resrart file for '//finp(1:mlen(finp)) call wrda (1, 1, kd, len(str), str) kd = iword * 200 + 100 kd_sta = kd kd_xy = kd call wrda (1, iword, kd, nxp, xm) call wrda (1, iword, kd, nyp, ym) kd_z = kd do i = 1, nz zz(i) = -zin(i) enddo call wrda (1, iword, kd, nz, zz) call segm_from_iox (npt, nseg, iox, iseg) kd_seg = kd call wrda (1, iword, kd, 2*nseg, iseg) kd_uv = kd nstart = 1 nlaststart = 1 nlast = nsteps nscpu = 0 nswll = 0 inf(1) = 120 inf(2) = nxp inf(3) = nyp inf(4) = nz inf(5) = 1 inf(6) = npt inf(7) = nseg inf(8) = ntrac inf(11) = nstart inf(15) = nlaststart inf(16) = nsteps inf(17) = nlast rnf(1) = delt rnf(2) = enso_start rnf(3) = enso_scale rnf(4) = rnf(2) return end subroutine read_rstrt (nxp, nyp, nz, npt) c-------------------------------------------------------------------- include 'comm_new.h' include 'comm_data.h' common /run/ nstart,nlaststart,nskip,nsteps,nergy,nskipo,nlast common /new_save/ iword, nruns, nscpu, nswll, inf(100), rnf(100) equivalence (inf(20),kd_sta), * (inf(21),kd_xy), (inf(22),kd_z), (inf(23),kd_time), * (inf(24),kd_seg), (inf(25),kd_uv), (inf(26),kd_temp), * (inf(27),kd_salt),(inf(28),kd_tr), (inf(29),kd_phi), * (inf(30),kd_conv), (inf(31),kd_means), * (inf(100),kd_end) iword = MACHINE_WORD call opda (1, 0, fbi(1:n_in)) kd = 0 call rdda (1, iword, kd, 100, inf) kd = irdda (iword, 100, rnf) ntype = idig(inf(1),1) if (ntype .ne. ITYPE_RST) then write (6, *) 'The file <' ,fbi(1:n_in),'> is not a RESTART-TYPE one !!' stop endif if (nxp.ne.inf(2) .or. nyp.ne.inf(3) .or. * nz.ne.inf(4) .or. npt.ne.inf(6) ) then write(6,*)'The model dimensions are different from the restart data!' stop endif nptz = inf(4)*inf(6) mtra = inf(8) nstart = inf(11) nruns = inf(12) nscpu = inf(13) nswll = inf(14) if (irest .eq. 1) then nlaststart = inf(15) else nlaststart = nstart endif nlast = nlaststart + nsteps - 1 enso_start = rnf(2) rnf(4) = rnf(2) + rnf(3)*nstart if (irest .ne. 3) then ekf1 = rnf(10) epf1 = rnf(11) hcf1 = rnf(12) wcf1 = rnf(13) vlf1 = rnf(14) do k = 1, inf(4)+1 hsave(k) = rnf(20+k) enddo endif kd = kd_uv call rdda (1, iword, kd, nptz, u) kd = irdda (iword, nptz, v) kd = irdda (iword, nptz, uc) kd = irdda (iword, nptz, vc) kd = irdda (iword, nptz, h) if (kd_temp .ne. 0) then kd = kd_temp call rdda (1, iword, kd, nptz, t) endif if (use_salt .and. kd_salt .ne. 0) then kd = kd_salt call rdda (1, iword, kd, nptz, sal) kd = irdda (iword, nptz, dens) endif if (use_trac .and. kd_tr .ne. 0) then kd = kd_tr call rdda (1, iword, kd, nptz*mtra, tr) endif if (kd_conv .ne. 0) then kd = kd_conv call rdda (1, iword, kd, nptz, convn) endif if (kd_means .ne. 0) then kd = kd_means call rdda (1, iword, kd, nptz, um) kd = irdda (iword, nptz, vm) kd = irdda (iword, nptz, tm) if (use_salt .and. kd_salt .ne. 0) then kd = irdda (iword, nptz, salm) kd = irdda (iword, nptz, densm) endif if (use_trac .and. kd_tr .ne. 0) kd = irdda (iword, nptz*mtra, trm) endif call clda(1) call opda (1, 1, fbo(1:n_out)) inf(8) = ntrac inf(12) = nruns + 1 inf(15) = nlaststart inf(16) = nsteps inf(17) = nlast rnf(1) = delt return end subroutine keep_rstrt (nstep, nskip) c-------------------------------------------------------------------- include 'comm_new.h' include 'comm_data.h' common /new_save/ iword, nruns, nscpu, nswll, inf(100), rnf(100) equivalence (inf(20),kd_sta), * (inf(21),kd_xy), (inf(22),kd_z), (inf(23),kd_time), * (inf(24),kd_seg), (inf(25),kd_uv), (inf(26),kd_temp), * (inf(27),kd_salt),(inf(28),kd_tr), (inf(29),kd_phi), * (inf(30),kd_conv), (inf(31),kd_means), * (inf(100),kd_end) if (nskip.eq.0) return if (mod(nstep, nskip) .ne. 0) goto 100 nptz = inf(4)*inf(6) kd = kd_uv call wrda (1, iword, kd, nptz, u) kd = iwrda (iword, nptz, v) kd = iwrda (iword, nptz, uc) kd = iwrda (iword, nptz, vc) kd = iwrda (iword, nptz, h) if ( use_temp ) then kd_temp = kd kd = iwrda (iword, nptz, t) endif if ( use_salt ) then kd_salt = kd kd = iwrda (iword, nptz, sal) kd = iwrda (iword, nptz, dens) endif if ( use_trac ) then kd_tr = kd kd = iwrda (iword, nptz*ntrac, tr) endif kd_conv = kd kd = iwrda(iword, nptz, convn) if ( save_mean ) then kd_means = kd kd = iwrda (iword, nptz, um) kd = iwrda (iword, nptz, vm) kd = iwrda (iword, nptz, tm) if ( use_salt ) then kd = iwrda (iword, nptz, salm) kd = iwrda (iword, nptz, densm) endif if ( use_trac ) kd = iwrda (iword, nptz*ntrac, trm) endif kd_end = kd inf(11) = nstep rnf(10) = ekf1 rnf(11) = epf1 rnf(12) = hcf1 rnf(13) = wcf1 rnf(14) = vlf1 do k = 1, inf(4)+1 rnf(20+k) = hsave(k) enddo 100 continue inf(18) = nstep inf(13) = nscpu + ipast_scpu() inf(14) = nswll + ipast_swll() kd = 0 call wrda (1, iword, kd, 100, inf) kd = iwrda (iword, 100, rnf) call flda (1) return end subroutine close_rstrt c--------------------------------- call clda(1) return end subroutine dump_rstrt c------------------------------------------------------------------ include 'comm_new.h' include 'comm_data.h' character*100 str character*4 mname(12) common /new_save/ iword, nruns, nscpu, nswll, inf(100), rnf(100) save mname data mname /'Jan.','Feb.','Mar.','Apr.','May ','Jun.', * 'Jul.','Aug.','Sep.','Oct.','Nov.','Dec.'/ iword = MACHINE_WORD call opda (1, 0, fbi(1:n_in)) kd = 0 call rdda (1, iword, kd, 100, inf) kd = irdda (iword, 100, rnf) kd = irdda (1, 100, str) call clda(1) write (6, *) 'Dump of the data file <',fbi(1:n_in),'>:' write (6, *) 'LABEL: ', str write (6, *) 'TYPE', inf(1) write (6, *) 'NX =', inf(2) write (6, *) 'NY =', inf(3) write (6, *) 'NZ =', inf(4) write (6, *) 'NT =', inf(5) write (6, *) 'NPACK =', inf(6) write (6, *) 'NSEGM =', inf(7) write (6, *) 'NTRACERS =', inf(8) call enso2date (rnf(2), id, im, iy) write (6, 101) 'Model Starting Date: ', mname(im), id, iy call enso2date (real(rnf(2)+rnf(3)*inf(17)), id, im, iy) write (6, 101) 'Scheduled End of the Run: ', mname(im), id, iy call enso2date (real(rnf(2)+rnf(3)*inf(18)), id, im, iy) write (6, 101) 'Current End of the Run: ' , mname(im), id, iy call enso2res (real(rnf(3)*inf(17)), id, im, iy) write (6, 201) 'Scheduled Length of the Run:', iy, im, id, inf(17) call enso2res (real(rnf(3)*inf(18)), id, im, iy) write (6, 201) 'Elapsed Model Time: ', iy, im, id, inf(18) write (6, *) 'Number of restarts:', inf(12) call enso2date (rnf(4), id, im, iy) write (6, 101) 'Last Re-start: ', mname(im), id, iy call enso2date (real(rnf(2) + rnf(3)*inf(15)), id, im, iy) write (6, 101) 'Last New Run: ', mname(im), id, iy call enso2date (real(rnf(2) + rnf(3)*inf(11)), id, im, iy) write (6, 101) 'Last Save: ', mname(im), id, iy i = inf(13) if (i/3600 .eq. 1) then write (6, 301) 'Total CPU time:', i/3600, mod(i,3600)/60, mod(i,60) else write (6, 302) 'Total CPU time:', i/3600, mod(i,3600)/60, mod(i,60) endif i = inf(14) if (i/3600 .eq. 1) then write (6, 301) 'Total WALL time:', i/3600, mod(i,3600)/60, mod(i,60) else write (6, 302) 'Total WALL time:', i/3600, mod(i,3600)/60, mod(i,60) endif 101 format (a30, a4, i2, ',' , i5) 201 format (a30, i5, ' Years', i3, ' months', i3, ' days. (', i7, ' steps.)') 301 format (a20, i5, ' hour ', i3, ' min', i3, ' sec.') 302 format (a20, i5, ' hours', i3, ' min', i3, ' sec.') return end subroutine segm_from_iox (npt, nseg, iox, iseg) c------------------------------------------------------ dimension iox(1), iseg(2,1) ista = iox(1) inext = ista + 1 nseg = 0 icount = 0 do i = 2, npt icurr = iox(i) icount = icount + 1 if (inext .ne. icurr) then nseg = nseg + 1 iseg(1,nseg) = ista iseg(2,nseg) = ista + icount - 1 ista = icurr icount = 0 endif inext = icurr + 1 enddo nseg = nseg + 1 iseg(1,nseg) = ista iseg(2,nseg) = ista + icount return end subroutine model_input(npt) c------------------------------------------------------------ c.....read model parameters and grid from file. implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' include 'comm_new.h' include 'comm_data.h' dimension dzin(MAXNZ+1) character*50 ubeg,vbeg,hbeg,tbeg,uout,vout,hout,tout,eout, + blk,wndx,wndy,cwndx,cwndy,cloud,ccloud,wout common/files/ubeg,vbeg,hbeg,tbeg,uout,vout,hout,tout,eout,wout common/run/ nstart,nlaststart,nskip,nsteps,nergy,nskipo,nlast common/param0/iyear,iday,isec,delt,ncyc,mbc,nonlin,label(20), + itherm,mlc,limp common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch common/winds/mtau,matau,tausc,atau,froude,icloud common/wnfils/wndx,wndy,cwndx,cwndy,cloud,ccloud common/shapi/nord,nshapu,nshaph,nshapt,mshapu,mshaph,mshapt common /vert/ zin(MAXNZ+1), hin(MAXNZ), t_in(MAXNZ+1), s_in(MAXNZ+1), * bint(MAXNZ), cint(MAXNZ), dzin(MAXNZ+1), sigma(MAXNZ) common/strech/xs(MAXXS),alpha(MAXXS),beta(MAXXS) common /errors/ ioerr, nstep common /new_filt/ MAXFO, nbxk, nbyk, nfxk, nfpxk, nfyk c------------------------------------------------------------------------------- common /pbl_param/ pbl_pnu,pbl_delta,pbl_pml,pbl_depth,pbl_betav, * pbl_grad,nstep_pbl, ipbl_advec, ipbl_jsta,ipbl_jend, * pbl_south,pbl_north, pbl_wmin c------------------------------------------------------------------------------- character*80 f_bar common /baro_files/ n_bar, f_bar integer iflag(12), RM12_NN real aflag(12) common /y12m_input/ RM12_NN, aflag, iflag common /baro_input/ n_def_cor, mod_scheme, mod_solver, BAR_DELTA, * BAR_DSINK, ibar_key, nbaro, rayl c------------------------------------------------------------------------------- real inp_flt, inp_days logical inp_def dimension flt(100) save flt character*80 fbdir external time, ctime integer time character*24 ctime c------------------------------------------------------------------------------- call inp_file(finp) call inp_vrnt(c_str(Variant), inp_int(c_str(Use_variant), -1)) if (inp_int(c_str(Trace),0) .ne. 0) * call inp_trace(finp(1:mlen(finp))//'.tr\0') n_dir= inp_str(c_str(Output_dir),c_str(output), fbdir) n_tios= inp_str(c_str(Base_file), * fbdir(1:n_dir)//'/'//finp(1:mlen(finp))//'\0',fbt) n_out = inp_str(c_str(Save_file), fbt(1:n_tios)//'.save\0', fbo) n_in = inp_str(c_str(Restart_file), fbo(1:n_out)//'\0', fbi) n_wnd = inp_str(c_str(Wind_file), c_str(wind_data), fbwnd) n_tem = inp_str(c_str(Temp_file) , c_str(temp_data), fbtem) n_sal = inp_str(c_str(Salt_file), c_str(salt_data), fbsal) n_sst = inp_str(c_str(SST_file), c_str(sst_data), fbsst) n_sss = inp_str(c_str(SSS_file), c_str(sss_data), fbsss) n_cld = inp_str(c_str(Cloud_file), c_str(cld_data), fbcld) n_slr = inp_str(c_str(Solar_file), c_str(slr_data), fbslr) n_evp = inp_str(c_str(EP_file), c_str(ep_data), fbevp) n_hcl = inp_str(c_str(MxlTcl_file),c_str(mltc_data), fbhcl) n_wsp = inp_str(c_str(Wndspd_file),c_str(wndspd_data), fwsp) n_uwd = inp_str(c_str(Uwnd_file),c_str(uwnd_data), fuwd) n_vwd = inp_str(c_str(Vwnd_file),c_str(vwnd_data), fvwd) n_ah = inp_str(c_str(Airhum_file),c_str(ahum_data), fah) n_at = inp_str(c_str(Airtem_file),c_str(atem_data), fat) n_dep = inp_str(c_str(Bath_file), c_str(bath_data), fbdep) if ( inp_def(c_str(Map_file)) ) * n_map = inp_str(c_str(Map_file), c_str(map_data), fbmap) ipre = inp_int(c_str(Pre_proc), 0) ! 0- normal run; 1- pre-processing irest = inp_int(c_str(Restart), 0) ! 0 - start; 1 - restart lev_err = inp_int(c_str(Error_level), 1) if (irest .eq. 0) then open (unit = iout, file = fout) if (lev_err .ge. 1) then write (iout, *) 'Run with control file <',finp(1:mlen(finp)),'>' write (iout, *) 'STARTED: ', ctime(time()) endif else open (unit = iout, file = fout, access='APPEND') if (lev_err .ge. 1) then write (iout, *) 'RE-STARTED: (code=',irest,' at:', ctime(time()) endif endif ncyc = inp_int(c_str(Ncycles), 4) delt = inp_days (c_str(Time_step), 1./24.) stpd = 1./delt minseg = 4 mtmp = stpd - ncyc if (mtmp .lt. 0) mtmp = 0 nsteps = nint(stpd * inp_days(c_str(Run_time), 365.)) + mtmp nskip = nint(stpd * inp_days(c_str(Save_step), 5.)) save_mean = inp_int(c_str(Save_mean), 1) .eq. 1 nergy = nint(stpd * inp_days(c_str(Energy_step), 1.)) call inp_date (c_str(Starting_date), 1,1,2000, imon,iday,iyear) enso_start = date2enso (iday, imon, iyear) enso_scale = delt * 12./365. land = 0 iglob = inp_int(c_str(Periodic), 0) mbc = inp_int(c_str(BC_type), 1) NXP = inp_int(c_str(NX), 10) NYP = inp_int(c_str(NY), 10) NZ = inp_int(c_str(NZ), 2) nsig= inp_int(c_str(Nsigma), 0) if (nsig.eq.1.or.nsig.eq.2.or.nsig.gt.nz.or.nsig.lt.0)then print*,'number of sigma layers must be 0 or between 3 and nz' stop endif if(max0(nxp,nyp).gt.MAXSID) call wspace('MAXSID', max0(nxp,nyp)) if(NZ .gt. MAXNZ) call wspace('MAXNZ', NZ) alat = inp_flt(c_str(South), 0.) blat = inp_flt(c_str(North), real(nyp-1)) alon = inp_flt(c_str(West), 0.) blon = inp_flt(c_str(East), real(nxp-1)) mgrid = inp_int(c_str(Grid_type), 1) nystrch = inp_int(c_str(Y_stretch_type), 1) if (nystrch.eq.1) then nsx = inp_rarr(c_str(X_stretch_pnts), 0, xs, xs) nsx = inp_rarr(c_str(X_alpha), 0, alpha, alpha) nsx = inp_rarr(c_str(X_beta), 0, beta, beta) nsy = inp_rarr(c_str(Y_stretch_pnts), 0, xs(nsx+1), xs(nsx+1)) nsy = inp_rarr(c_str(Y_alpha), 0, alpha(nsx+1), alpha(nsx+1)) nsy = inp_rarr(c_str(Y_beta), 0, beta(nsx+1), beta(nsx+1)) if(nsx+nsy.gt.maxxs) call wspace('MAXXS', nsx+nsy) endif ibaro = inp_int(c_str(Baro_solv), 0) if (ibaro .ne. 0) then nbaro = inp_int(c_str(Baro_step), 6) n_def_cor = inp_int(c_str(Baro_defcor_step), 0) mod_scheme = inp_int(c_str(Baro_scheme), 9) mod_solver = inp_int(c_str(Baro_solver), 1) BAR_DSINK = inp_flt(c_str(Baro_depsink), 500.) BAR_DELTA = inp_flt(c_str(Baro_delta), .001) rayl = inp_flt(c_str(Baro_rayl), .002) ibar_key = inp_int(c_str(Baro_err), 1) if (ibar_key .ne. 0) * n_bar = inp_str(c_str(Baro_errfile), * finp(1:mlen(finp))//'.bar\0', f_bar) c------------FROM ".y12m"-------------------------------------------- c RM12_NN : size of workspace in numbers of NONZ c iflag(11) : number of iterations c iflag : see y12m documentation c aflag : see y12m documentation c default values: do i = 1, 12 iflag(i) = 0 aflag(i) = 0 enddo RM12_NN = inp_int(c_str(Baro_nnonz), 4) iflag(11) = inp_int(c_str(Baro_niter), 100) iflag(1) = 1 iflag(2) = 3 iflag(3) = 1 iflag(5) = 2 aflag(1) = 2. aflag(2) = 1.e-4 aflag(3) = 1.e6 aflag(4) = 1.e-12 aflag(5) = -1000. call inp_iarr (c_str(Baro_iflag), 5, iflag, iflag) call inp_rarr (c_str(Baro_aflag), 5, aflag, aflag) endif nonlin= inp_int(c_str(Nonlin), 1) if (nonlin.ne.1) then print*, 'Nonlin disabled (ie, always nonlinear), please do not use' stop endif itemp = inp_int(c_str(Thermo), 1) if (itemp.eq.0) then print*, 'Thermo = 0 disabled' stop endif isalt = inp_int(c_str(Salinity), 0) isolrp = inp_int(c_str(Solar_penrad), 0) if (isolrp .eq. 1) solr_gamma = inp_flt(c_str(Solar_gamma), .333333) icl_h = inp_int(c_str(Clim_H), 0) ! 0-init from hin;1-Cnst;2-Vary; if (icl_h .ne. 0.and.nsig.eq.0) then print*,'must use Nsigma > 2 if Clim_H <> 0' stop endif icl_htop = inp_int(c_str(Clim_htop), 0) ! 0-Const; 1-Vary; icl_ts = inp_int(c_str(Clim_TS), 0) ! 0-init from Tin;1-Cnst;2-Vary; if (icl_h.eq.2 .and. icl_ts.ne.2) * write (ioerr, *) 'Warning: check Clim_H vs. Clim_TS...' if (icl_ts.eq.0 .and. initt.ne.0) * write (ioerr, *) 'Warning: to initialize temp, use icl_ts <> 0' if (icl_ts.eq.0 .and. inits.ne.0) * write (ioerr, *) 'Warning: to initialize salt, use icl_ts <> 0' icl_rlx = inp_int(c_str(Clim_relax), 0) ! 0-no relax; 1-relax N-S; ksponge = inp_flt(c_str(Sponge_width), 5) clm_coef = inp_flt(c_str(Clim_coef), 0.1) ! Relaxation coefficient clm_no = inp_flt(c_str(Clim_nlat), 90.) ! North relaxation latitude. clm_so = inp_flt(c_str(Clim_slat), -90.) ! South relaxation latitude. imix = inp_int(c_str(Mixing), 0) if (imix .ne. 0) limp = inp_int(c_str(Mix_step), 3*ncyc) cm_mix = inp_flt(c_str(Mix_cm), 1.25) ! "m" turbulence coefficient cn_mix = inp_flt(c_str(Mix_cn), 0.17) ! "n" turbulence coefficient cn_mix = 1. - cn_mix hmin_mix = inp_flt(c_str(Mix_hmin), 10.) ! H1 - min hmax_mix = inp_flt(c_str(Mix_hmax), 100.) ! H1 - max ric1_mix = inp_flt(c_str(Mix_ric1), .65) ! Ri for k = 1 ric2_mix = inp_flt(c_str(Mix_ric2), .25) ! Ri for k = 2:NZ iuse_gam = inp_int(c_str(Mix_usegam), 0) ! 1 - use gammas in jpmix gam1_mix = inp_flt(c_str(Mix_gam1), 1.) ! Gamma for k = 1 gam2_mix = inp_flt(c_str(Mix_gam2), 1.) ! Gamma for k = 2:NZ iwnd_mix = inp_int(c_str(Mix_wnd), 0) ! 0:use tau; 1: use windspeeds if (initq .ne. 8) iwnd_mix = 0 ! valid only if PBL is "ON" mix_wtop = inp_int(c_str(Mix_wtop), 2) if (imix .eq. 4.and.nsig .eq. 0) then print*, 'must set Nsigma>1 to use variable depth mixed layer' stop endif ntrac = inp_int(c_str(Tracers), 0) use_temp = (itemp .ne. 0) use_salt = (isalt .ne. 0) use_trac = (ntrac .ne. 0) iv_top = inp_int(c_str(Vert_Top), 1) if (iv_top.ne.1.and.nsig.eq.0) then print*, 'must use Sigma layers for non-constant depth mixed layer' endif iv_bot = inp_int(c_str(Vert_Bot), 99) !! 1-CNST;2-Sigm;3-W=0;4-No Motion if (iv_bot.ne.99) then print*, 'Vert_Bot ignored, use Nsigma to determine last layer type' endif iv_bump = inp_int(c_str(Vert_Bump), 0) !! 0-Sigm;2-Sm.bmp;3-Lg.bmp;4-Exp. if (iv_bump.ne.0) then print*, 'no interface bumps allowed, reset Vert_Bump = 0' stop endif iv_fix = inp_int(c_str(Vert_Fix), 1) !! 0-Free Surface; 1-Fixed Depth if (iv_fix.ne.1) then print*, 'reduced gravity runs disabled, reset Vert_Fix = 1' stop endif iv_sys = inp_int(c_str(Vert_Sys), 99) !! disabled if (iv_sys.ne.99) then print*, 'Vert_Sys ignored, H is filtered only in sigma layers' endif mbot_bc = inp_int(c_str(BC_bot_temp), 0) !! 0-constant temp; 1 - d/dz=0 if (mbot_bc.ne.1) then print*, 'do you REALLY want constant bottom temp with bathymetry?' print*, 'reconsider using BC_bot_temp = 1' stop endif nordu = inp_int(c_str(Shap_vel_order), 4) mshapu = inp_int(c_str(Shap_vel_type), 3) nshapu = inp_int(c_str(Shap_vel_step), 3*ncyc) nordh = inp_int(c_str(Shap_scl_order), 4) mshaph = inp_int(c_str(Shap_scl_type), 1) nshaph = inp_int(c_str(Shap_scl_step), 3*ncyc) nord = max(nordu,nordh) MAXFO = nord mtau = inp_int(c_str(Wind_forc), 0) if (mtau .eq. 5) itau_cos = inp_int(c_str(Wind_cos), 0) call inp_rarr (c_str(Wind_tauxy), 2, flt, flt) tausc = flt(1) atau = flt(2) initt = inp_int(c_str(Temp_init), 0) if (initt.ne.0.and.initt.ne.3) then print*,'must use Temp_init = 0 or 3' stop endif inits = inp_int(c_str(Salt_init), 0) if (inits.ne.0.and.inits.ne.3) then print*,'must use Salt_init = 0 or 3' stop endif froude = inp_flt(c_str(Temp_froude), 0.01) initq = inp_int(c_str(Heat_forc), 0) initep = inp_int(c_str(EP_forc), 0) icloud = inp_int(c_str(Cloud_forc), 0) qcon = inp_flt (c_str(Rho_CP), 4.12e6) rlx_time = 86400.*inp_days (c_str(Rlx_time), 30.) if (initq .ge. 8) then nstep_pbl = nint(stpd * inp_days(c_str(PBL_step), 5.)) ipbl_advec = inp_int(c_str(PBL_advec), 1) pbl_wmin = inp_flt(c_str(PBL_wmin), 4.) pbl_pnu = inp_flt(c_str(PBL_pnu), 0.4e+7) pbl_delta = inp_flt(c_str(PBL_delta), 0.25 ) pbl_pml = inp_flt(c_str(PBL_pml), 6000.) pbl_depth = inp_flt(c_str(PBL_depth), 600. ) pbl_betav = inp_flt(c_str(PBL_betav), 0.17 ) pbl_grad = inp_flt(c_str(PBL_grad), -2.) pbl_south = inp_flt(c_str(PBL_south), -200.) pbl_north = inp_flt(c_str(PBL_north), 200.) endif iherr = inp_int(c_str(Layer_volume), 0) initb = inp_int(c_str(Bathymetry), 0) initbt = inp_int(c_str(Bath_type), 0) dep_min = inp_flt(c_str(Bath_min), 100) if (nz .gt. inp_rarr(c_str(Z_profile), nz, zin, zin).and. * nz .gt. inp_rarr(c_str(H_profile), nz, hin, hin)) then call perror1('Z_profile or H_profile: not enough terms...',-1) stop endif zin(1) = hin(1)/2. z_bot = hin(1) do k = 2, nz zin(k) = - zin(k-1) + 2.*z_bot z_bot = z_bot + hin(k) enddo zin(nz+1) = z_bot dzin(1) = zin(1) do j = 2, nz dzin(j) = (zin(j)-zin(j-1))/2. enddo dzin(nz+1) = zin(nz+1)-zin(nz) c Z_profile specification overrides H_profile zdum = inp_flt(c_str(Z_bot), 0) if (zdum.eq.0.and.z_bot.eq.0) then print*,'must specifiy Z_bot as well as Z_profile' stop elseif (zdum.gt.0) then c reread Z_profile ndum = inp_rarr(c_str(Z_profile), nz, zin, zin) zin(1) = zin(2)/3. ! enforce this condition zin(nz+1) = zdum dzin(1) = zin(1) do j = 2, nz dzin(j) = (zin(j)-zin(j-1))/2. if (dzin(j).lt.0.)then print*,'negative layer depths, check Z_profile' stop endif enddo dzin(nz+1) = zin(nz+1)-zin(nz) if (dzin(nz+1).lt.0.)then print*,'negative depth of last layer, check Z_bot' stop endif hin(1) = 2.*dzin(1) do j = 1, nz hin(j) = dzin(j) + dzin(j+1) enddo endif if (hin(1).eq.0) then print*,'must specify either Z_profile or H_profile' stop endif z_begin = zin(nsig) sigma(1) = -1./3. sigma(2) = sigma(1) do j = 3, nsig sigma(j) = dzin(j)/(z_begin-3.*dzin(1)) enddo do j = nsig + 1, nz sigma(j) = 0. enddo trans_coef = (zin(1)+zin(2))/2./rlx_time i = inp_rarr(c_str(T_profile), nz, t_in, t_in) if (i .lt. nz) then call perror1('T_profile: not enough terms...',-1) elseif (i .eq. nz) then t_in(nz+1) = inp_flt(c_str(Temp_bot), 0.) endif TATM = inp_flt(c_str(Temp_atm), 30.) i = inp_rarr(c_str(S_profile), nz, s_in, s_in) if (i .lt. nz) then call perror1('S_profile: not enough terms...',-1) elseif (i .eq. nz) then s_in(nz+1) = inp_flt(c_str(Salt_bot), 36.) endif SATM = inp_flt(c_str(Salt_atm), 35.4) TEMP_BOT = t_in(nz+1) SALT_BOT = s_in(nz+1) POTND_BOT = inp_flt(c_str(Dens_bot), pdens_pnt (TEMP_BOT, SALT_BOT)) c.....Bi & Ci are now *real* "Nu" & "Ka" (not scaled by depth): bi = inp_flt (c_str(Bint), 1.e-3) ci = inp_flt (c_str(Cint), 1.e-4) do i = 1, nz bint(i) = bi cint(i) = ci enddo call inp_str(c_str(Grid_label),'Data: '//finp(1:mlen(finp))//'\0',label) call mem_alloc (p_mask, NXP*NYP*NZ, 1, 'mask') if (n_map .ne. 0) call inp_file(fbmap(1:n_map)//'\0') call read_mask ('Grid_mask', nxp, nyp, nxyc, mask) call mem_alloc (p_iox, nxyc, 1, 'iox') npt = nxyc npt1 = 1 npt2 = 1 + npt npt3 = 1 + 2*npt npt4 = 1 + 3*npt first_step = .true. if (lev_err .ge. 1 .and. irest .eq. 0) then write (iout, *) 'NX =', nxp write (iout, *) 'NY =', nyp write (iout, *) 'NZ =', nz write (iout, *) 'NPACK =', nxyc write (iout, *) 'NTRACERS =', ntrac if (ibaro .eq. 0) then write (iout, *) 'BARO OFF' else write (iout, *) 'BARO ON' endif call flush(iout) endif return end subroutine read_mask (tag, nx, ny, npack, mask) c----------------------------------------------------- character*(*) tag dimension mask(nx, 1) character*1 ch, buff(119) character*64 number logical inp_def equivalence (ch, buff(1)), (number, buff(2)) if ( inp_def(tag//'\0') ) then ix = 1 jy = ny ic0 = ichar('0') do while ( inp_wnxt(buff) .gt. 0) if (ch.eq.'0' .or. ch.eq.'1' .or. ch.eq.'2') then do i = 1, nx mask(i,jy) = ichar(buff(i)) - ic0 enddo jy = jy - 1 elseif (ch.eq.'w'.or.ch.eq.'z'.or.ch.eq.'x'.or.ch.eq.'s') then read (number, *) kx ix = ix + kx if (ix - 1 .gt. NX) goto 200 if (ch.eq.'w') ich = 1 if (ch.eq.'x') ich = 2 if (ch.eq.'s') ich = 3 if (ch.eq.'z') ich = 0 do i = ix-kx, ix-1 mask(i,jy) = ich enddo if (ix-1 .eq. NX) then ix = 1 jy = jy - 1 endif elseif (ch .eq. 'r') then read (number, *) ky jy = jy - ky + 2 if (jy .lt. 1) goto 200 do j = jy+ky-2, jy, -1 do i = 1, nx mask(i,j) = mask(i,j+1) enddo enddo jy = jy - 1 else goto 200 endif if (jy .eq. 0) goto 100 enddo 100 npack = 0 do j = 1, ny do i = 1, nx if (mask(i,j) .ne. 0) npack = npack + 1 enddo enddo else npack = nx*ny do j = 1, ny do i = 1, nx mask(i,j) = 1 enddo enddo endif return 200 write (6, *) '!!!read_mask: wrong mask data, i,j:', ix,jy stop end