#define ITYPE_RST 0 #define MACHINE_WORD 4 #define c_str(s) ('s\0') #define n_a(s,i) ('s'//numbers(i:i)//'\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_mask),(inf(32),kd_psi), (inf(33),kd_ice), * (inf(100),kd_end) call opda (1, 1, fbo(1:n_out)) call opda (0, 1, fbo0(1:n_out0)) do i = 1, 100 inf(i) = 0 rnf(i) = 0. enddo iword = MACHINE_WORD kd = iword * 200 str = 'Restart file for '//finp(1:mlen(finp)) call wrda (1, 1, kd, len(str), str) call wrda (0, 1, kd, len(str), str) kd = iword * 200 + 100 kd_sta = kd kd_xy = kd kd_keep = kd_xy call wrda (1, iword, kd, nxp, xm) call wrda (1, iword, kd, nyp, ym) kd_xy = kd_keep call wrda (0, iword, kd, nxp, xm) call wrda (0, iword, kd, nyp, ym) kd_z = kd kd_keep = kd_z do i = 1, nz zz(i) = -zin(i) enddo call wrda (1, iword, kd, nz, zz) kd_z = kd_keep call wrda (0, iword, kd, nz, zz) call segm_from_iox (npt, nseg, iox, iseg) kd_seg = kd kd_keep = kd_seg call wrda (1, iword, kd, 2*nseg, iseg) kd_seg = kd_keep call wrda (0, 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_tot 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) budget_sp = 0. budget_ri = 0. budget_re = 0. budget_eb = 0. budget_pr = 0. coef_bud = coef_precip return end subroutine read_rstrt (nxp, nyp, nz, npt) c-------------------------------------------------------------------- include 'comm_new.h' include 'comm_data.h' include 'comm_bio.h' include 'diffiso.h' common /run/ nstart,nlaststart,nskip,nsteps,nergy,nskipo,nlast common /new_save/ iword, nruns, nscpu, nswll, inf(100), rnf(100) common /mean_comm/ nmcount, nmcount_bio, nmcount_mosf, nmcount_conv character*4 mname(12) save mname 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_mask),(inf(32),kd_psi), (inf(33),kd_ice), * (inf(100),kd_end) 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) 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) nptzp = inf(4)*inf(6) +inf(6) nruns = inf(12) nmcount = 0 nmcount_bio = 0 nmcount_mosf = 0 nmcount_conv = 0 budget_sp = 0. budget_ri = 0. budget_eb = 0. budget_pr = 0. coef_bud = coef_precip budget_re = 0. if (irest.eq.1) coef_bud = rnf(19) f0 = rnf(5) if (irest.ne.4) then nstart = inf(11) + 1 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 else nstart = 1 nlaststart = 1 nlast = nsteps nscpu = 0 nswll = 0 inf(11) = nstart inf(15) = nlaststart inf(16) = nsteps inf(17) = nlast rnf(2) = enso_start rnf(4) = rnf(2) endif if (irest .lt. 3 .and. irest. ge. 0) 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 (ibaro.ne.0.and.kd_psi .ne. 0) then kd = kd_psi call rdda (1, iword, kd, npt, psi) endif if (use_ice .and. kd_ice .ne. 0) then kd = kd_ice call rdda (1, iword, kd, npt, trs(i_hice)) kd = irdda (iword, npt, trs(i_cice)) kd = irdda (iword, npt, trs(i_thice)) endif if (kd_salt .ne. 0) then kd = kd_salt call rdda (1, iword, kd, nptz, sal) kd = irdda (iword, nptz, dens) endif if ((use_trac .or.use_bio) .and. kd_tr .ne. 0) then kd = kd_tr call rdda (1, iword, kd, nptz*mtra, tr) endif call clda(1) call opda (1, 1, fbo(1:n_out)) call opda (0, 1, fbo0(1:n_out0)) inf(8) = ntrac_tot inf(12) = nruns + 1 inf(15) = nlaststart inf(16) = nsteps inf(17) = nlast rnf(1) = delt call enso2date (rnf(2), id, im, iy) write (iout, 101) 'Model Starting Date: ', mname(im), id, iy call enso2date (real(rnf(2)+enso_scale*nlast), id, im, iy) write (iout, 101) 'Scheduled End of the Run: ', mname(im), id, iy call enso2date (real(rnf(2)+enso_scale*nstep), id, im, iy) write (iout, 101) 'Current Time: ' , mname(im), id, iy call enso2res (real(rnf(3)*inf(17)), id, im, iy) write (iout, 201) 'Scheduled Length of the Run:', iy, im, id, inf(17) call enso2res (real(rnf(3)*inf(18)), id, im, iy) write (iout, 201) 'Elapsed Model Time: ', iy, im, id, inf(18) write (iout, *) 'Number of restarts:', inf(12) 101 format (a30, a4, i2, ',' , i5) 201 format (a30, i5, ' Years', i3, ' months', i3, ' days. (', i7, ' steps.)') 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 c------------------------------------------------------------ subroutine model_input(npt,mbox) 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' include 'comm_bio.h' include 'biology.h' include 'amlice.h' include 'diffiso.h' include 'comm_topo.h' 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 common/wnfils/wndx,wndy,cwndx,cwndy,cloud,ccloud common /vert/ zin(MAXNZ+1), hin(MAXNZ), t_in(MAXNZ+1), s_in(MAXNZ+1), * bint(MAXNZ), cint(MAXNZ), dzin(MAXNZ+1), sigma(MAXNZ), * facz(MAXNZ), ttr_fac(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_DEPMIN, * BAR_DSINK, ibar_key, nbaro, rayl, nonlin_baro character*1 BC_I(0:9) common /masks/BC_I logical use_per_island, use_stan_island, use_sunk common /baro_island/ * alons_min(10),alons_max(10),alats_min(10),alats_max(10) * ,alon1_min,alon1_max,alat1_min,alat1_max * ,per_lat,use_sunk,use_per_island,use_stan_island * ,itrans,ftrans common /new_island/n_islands,in_rhs(0:9),ityp_island(0:9), * alon_island(9),alat_island(9),btrans(0:9), * n_poly(9),alon_poly(100,9),alat_poly(100,9), * i_poly(100,9),j_poly(100,9) c-------------------------------------------------------------- logical sink_region common /bathymetry/sink_region,sr_depth, alon_sr, blon_sr, * alat_sr, blat_sr logical use_hi common /order/ use_hi common/friction/ ibfric,b_fric, sw_fric 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 character*30 date,pwd character*10 name 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') call version (date,pwd,name) iret = inp_str(c_str(Compile_date), date//'\0', date) iret = inp_str(c_str(Compile_dir), pwd//'\0', pwd) iret = inp_str(c_str(Compile_who), name//'\0', name) n_debug = inp_int(c_str(Debug),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_out0 = inp_str(c_str(Save_file0), fbt(1:n_tios)//'.bak\0', fbo0) n_in = inp_str(c_str(Restart_file), fbo(1:n_out)//'\0', fbi) n_tau = inp_str(c_str(Wind_file), c_str(wind_data), fbtau) 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_psi = inp_str(c_str(Psi_file), c_str(psi_data), fbpsi) 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_prp = inp_str(c_str(Precip_file), c_str(prp_data), fprp) 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) use_wndspd = fwsp(1:n_wsp).ne.'NONE' 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) n_sp = inp_str(c_str(Sponge_file), c_str(sponge_data), fbsp) n_riv = inp_str(c_str(River_file), c_str(river_data), fbriv) n_grad = inp_str(c_str(Grad_file), c_str(grad_data), fbgrad) if ( inp_def(c_str(Map_file)) ) * n_map = inp_str(c_str(Map_file), c_str(map_data), fbmap) irest = inp_int(c_str(Restart), 0) ! 0 - start; 1 - restart lev_err = inp_int(c_str(Error_level), 1) ncyc = inp_int(c_str(Ncycles), 4) delt = inp_days (c_str(Time_step), 1./24.) t_fac = inp_flt(c_str(Time_fac), 1.) bt_fac = inp_flt(c_str(Time_fac_baro), t_fac) stpd = 1./delt tr_zmin = inp_flt(c_str(Time_trace_zmin), 10000.) tr_zmax = inp_flt(c_str(Time_trace_zmax), 0.) tr_fac = inp_flt(c_str(Time_trace_fac), 1.) 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.)) nskip = max(4,nskip - mod(nskip,4)) 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. 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') write (iout, *)'---------------------------------------------------' if (lev_err .ge. 1) then write (iout, *) 'RE-STARTED: (code=',irest,' at:', ctime(time()) endif endif land = 0 iglob = inp_int(c_str(Periodic), 0) iperx = (iglob.eq.1) mbc = inp_int(c_str(BC_type), 2) if (mbc.eq.1.and.ibaro.ne.0) then print*,'using no-slip b.c. and barotropic, retry with free-slip b.c.' endif bc_time = inp_days (c_str(BC_coef), 30.) bc_coef = 0. if (bc_time.ne.0.) bc_coef = delt/bc_time mseg = inp_int(c_str(Discretization), 2) use_hi = (mseg .eq. 4) 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) mbox = inp_int(c_str(Box_type), 1) 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) nystrch = inp_int(c_str(Y_stretch_type), 1) if (nystrch.eq.1) then 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 ipole = inp_int(c_str(Pole_shift), 0) if (ipole .eq. 1) then pole_alp = inp_flt(c_str(Pole_alpha), 0.) pole_bet = inp_flt(c_str(Pole_beta), -90.) pole_gam = inp_flt(c_str(Pole_gamma), 0.) call comp_rotma (pole_alp, pole_bet, pole_gam) endif b_fric = inp_flt(c_str(Bottom_friction), .001) sw_fric = inp_flt(c_str(Sidewall_friction), b_fric) ibfric = inp_int(c_str(Bottom_draglaw), 1) dep_min = inp_flt(c_str(Bath_min), 100.) ibaro = inp_int(c_str(Baro_solv), 0) if (ibaro .ne. 0) call baro_init 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 idens = inp_int(c_str(Density), 1) isalt = inp_int(c_str(Salinity), 4) 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; icl_psi = inp_int(c_str(Clim_psi),0) ! 0-don't use;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; nz_fix = inp_int(c_str(Clim_nz), nz+1) icl_fix = inp_int(c_str(Clim_nz_type), 1) 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.gt.1) * write (ioerr, *) 'Warning: to initialize temp, use icl_ts <> 0' if (icl_ts.eq.0 .and. inits.gt.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_int(c_str(Sponge_width), 5) krelax = inp_int(c_str(Relax_width), 5) clm_time = inp_days (c_str(Sponge_time), 30.) clm_time_v = inp_days (c_str(Sponge_time_v), 10000.) clm_time_tr = inp_days (c_str(Sponge_time_tr), 0.) clm_time_trb = inp_days (c_str(Sponge_time_trb), 0.) clm_time_psi = inp_days (c_str(Sponge_time_psi), 30.) clm_coef = ncyc*delt/clm_time clm_coef_v = clm_time/clm_time_v clm_coef_tr = 0. clm_coef_trb = 0. if (clm_time_tr.ne.0) clm_coef_tr = ncyc*delt/clm_time_tr if (clm_time_trb.ne.0) clm_coef_trb = ncyc*delt/clm_time_trb clm_psi = clm_coef* nbaro*clm_time/clm_time_psi 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.eq.8) then ! kpp mixing kpp_type = inp_int(c_str(KPP_type), 0) if (kpp_type.eq.0) print*,'Using windspeed version of KPP mixing' if (kpp_type.eq.1) print*,'Using windstress version of KPP mixing' endif 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 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) iice = inp_int(c_str(Ice), 0) mice = inp_int(c_str(Ice_dynamics), 0) if (iice.eq.0) mice = 0 use_trac = (ntrac .ne. 0) use_ice = (iice .ne. 0) use_dyice_old = (mice .eq. 1) ! crude ice dynamics use_dyice = (mice .eq. 2) ! Bruno's ice dynamics ibio_type = inp_int(c_str(Bio_type),0) use_bio = ibio_type .ne. 0 use_bio_old = ibio_type .eq. 2 use_river = (inp_int(c_str(Rivers), 0) .eq. 1) coef_river = inp_flt(c_str(River_coef),1.) dmin_river = inp_flt(c_str(River_dmin),5.) salt_fixed = inp_flt(c_str(Salt_fixed),36.) ntrac_bio = 0 if (use_bio) ntrac_bio = inp_int(c_str(Bio_ntrac),NTRBIOdef) ntrac_tot = ntrac + ntrac_bio ibio_key = inp_int(c_str(Bio_err), 0) if (ibio_key.ne.0) * n_bio = inp_str(c_str(Bio_errfile), * finp(1:mlen(finp))//'.bio\0', f_bio) nstep_bio = max(1,nint(stpd * inp_days(c_str(Bio_step), delt))) nstep_bio_print = nint(stpd * inp_days(c_str(Bio_print), 1.)) nstep_propag = nint(stpd * inp_days(c_str(Bio_step_propag), 1.)) if (ntrac_tot .gt. ntrac_max) then print*, 'need to increase NTRAC_MAX in comm_new.h' stop endif if (use_bio) call bio_init do i = 1, ntrac tr_in(i) = 100. enddo if (use_bio) then inittr(ntrac+ndet) = 2 inittr(ntrac+nzoo) = 2 tr_in(ntrac+ndet) = 0.1 tr_in(ntrac+nzoo) = 0.0143 inittr(ntrac+nnut) = 4 inittr(ntrac+nphy) = 4 endif nret = inp_iarr (c_str(Trac_init), ntrac_tot, inittr, inittr) nret = inp_rarr (c_str(Trac_init_value), ntrac_tot, tr_in, tr_in) if (irest.gt.0) then do i = 1, ntrac_tot inittr(i) = 0 enddo endif if (use_bio.or.use_trac) then do i = 1, ntrac ftrnm_def(i) = "TR" enddo ftrnm_def(ntrac+ndet) = "DET" ftrnm_def(ntrac+nzoo) = "ZOO" ftrnm_def(ntrac+nphy) = "PHY" ftrnm_def(ntrac+nnut) = "NUT" endif num1 = inp_sarr(c_str(Tr_name),ntrac_tot,ftrnm_def,80,name_tr,ftrnm) num2 = inp_sarr(c_str(Tr_clim_file),0,fbtr,80,n_tr,fbtr) if_tr_adv = inp_int(c_str(Tr_adv_type),1) if (if_tr_adv.eq.0) print*,'Tracer advection turned off' if (use_bio) then itra = ntrac + nnut n_tr(itra)=inp_str(c_str(Nitrate_file),c_str(nitrate_data),fbtr(itra)) endif c write(*,*) " " c write(*,*) " ftrnm= ", ftrnjm c write(*,*) " fbtr= ", fbtr c write(*,*) " tr_atm= ", tr_atm c write(*,*) " number1= ", num1 c write(*,*) " number2= ", num2 isod = inp_int(c_str(Iso_diffusion), 0) use_trdiff = (isod .ne. 0) isod = inp_int(c_str(Diffiso_scl), 0) use_trdiff = use_trdiff.or.(isod .ne. 0) isod = inp_int(c_str(Diffiso_vel), 0) use_modiff = (isod .ne. 0) use_diffiso = use_trdiff.or.use_modiff if (use_modiff.and.nsig.ne.0) then print*, 'cannot use non-shapiro diffusion options with sigma layers' endif 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), 1) !! 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) if (mshapu.lt.0) then print*,'Shap_vel_type reset to 0' mshapu = 0 endif nshapu = inp_int(c_str(Shap_vel_step), 3*ncyc) dshapu = inp_flt(c_str(Shap_vel_damp), 1.0) dshapv = inp_flt(c_str(Shap_vel_dampv), dshapu) nordh = inp_int(c_str(Shap_scl_order), 4) mshaph = inp_int(c_str(Shap_scl_type), 0) nshaph = inp_int(c_str(Shap_scl_step), 3*ncyc) dshaph = inp_flt(c_str(Shap_scl_damp), 1.0) nordp = inp_int(c_str(Shap_psi_order), 4) mshapp = inp_int(c_str(Shap_psi_type), 0) nshapp = inp_int(c_str(Shap_psi_step), 1) dshapp = inp_flt(c_str(Shap_psi_damp), 1.0) nord = max(nordu,nordh,nordp,nordd) MAXFO = nord mtau = inp_int(c_str(Wind_forc), 0) if (mtau .eq. 5) itau_cos = inp_int(c_str(Wind_cos), 0) flt(1) = 1. flt(2) = 1. nret = inp_rarr (c_str(Wind_tauxy), 2, flt, flt) tausc = flt(1) atau = flt(2) initt = inp_int(c_str(Temp_init), 0) temp_coef = inp_flt(c_str(Temp_coef), 0.) 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) if (initep.eq.8) then coef_precip = inp_flt(c_str(EP_coef_precip), 1.) v_sponge = inp_flt(c_str(EP_sponge), 0.) use_budget = inp_int(c_str(EP_budget), 0) .ne. 0 tscl_bud = inp_days (c_str(EP_bud_time), 0) if (tscl_bud.ne.0) tscl_bud = ncyc*delt/tscl_bud print*,'budget_coef (tscl_bud),coef ',tscl_bud,coef_precip endif iep_key = inp_int(c_str(EP_err), 0) if (iep_key .ne. 0) * n_ep = inp_str(c_str(EP_errfile), * finp(1:mlen(finp))//'.ep\0', f_ep) if (iep_key.gt.0) then if (irest .eq. 1. .or. irest .eq. 2) then open (unit = IEP_OUT, file = f_ep(1:n_ep), access='APPEND') else open (unit = IEP_OUT, file = f_ep(1:n_ep)) write(IEP_OUT,*)' Salt Budget' write(IEP_OUT,*) write(IEP_OUT,*)' coef_pr evap-brne precip sponge river relax total' write(IEP_OUT,*) endif endif igas = inp_int(c_str(Gas_exchange), 0) use_wnsp = (igas .ne. 0) if (use_ice.or.initq.eq.8) use_wnsp = .true. if (use_ice.and.initq.ne.8) then print*,'Using AML-ICE, so I am ignoring your Heat_forc parameter' endif if (initq.ne.8 .and. iice.ne.1)iwnd_mix = 0 ! valid only if PBL is "ON" qcon = inp_flt (c_str(Rho_CP), 4.12e6) tdef = 30. if (use_ice) tdef = 0. rlx_time_sst = 86400.*inp_days (c_str(Rlx_time_sst), tdef) rlx_time_sss = 86400.*inp_days (c_str(Rlx_time_sss), 30.) itank = inp_int(c_str(Tank_type), 0) if (itank.ne.0) cint_fac = inp_int(c_str(Tank_cint_fac), 0.) if (itank.ne.0) bint_fac = inp_int(c_str(Tank_bint_fac), 0.) if (initq .ge. 8 .or. use_ice) 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.) pbl_eddy = inp_flt(c_str(PBL_eddy), pbl_eddydef) pbl_eddy_nlat = inp_flt(c_str(PBL_eddy_nlat),pbl_eddy_nlat_def) pbl_eddy_nwidth= inp_flt(c_str(PBL_eddy_nwidth),pbl_eddy_nwidth_def) pbl_eddy_slat = inp_flt(c_str(PBL_eddy_slat),pbl_eddy_slat_def) pbl_eddy_swidth= inp_flt(c_str(PBL_eddy_swidth),pbl_eddy_swidth_def) pbl_eddy_xbox= inp_flt(c_str(PBL_eddy_xbox),pbl_eddy_xbox_def) pbl_eddy_ybox= inp_flt(c_str(PBL_eddy_ybox),pbl_eddy_ybox_def) c0_fac = inp_flt(c_str(PBL_c0_fac), c0_facdef) endif mshapd = 0 initkmap = inp_int(c_str(Diffiso_map), 0) if (use_diffiso) then c default values are given diffiso.h n_dif = inp_str(c_str(Dif_file), c_str(dif_data), fbdif) diffiso_alpha = inp_flt(c_str(Diffiso_alpha), diffiso_alphadef) if (diffiso_alpha.ne.0) then nordd = inp_int(c_str(Shap_dens_order), 4) mshapd = inp_int(c_str(Shap_dens_type), -1) nshapd = inp_int(c_str(Shap_dens_step), 1) dshapd = inp_flt(c_str(Shap_dens_damp), 1.0) endif diffiso_order = inp_int(c_str(Diffiso_order), 1) nslope = inp_int(c_str(Diffiso_step), 1) nkmap = 4*nint(stpd * inp_days(c_str(Kmap_step), 30.)/4.) diffiso_eps = inp_flt(c_str(Diffiso_eps), diffiso_epsdef) diffiso_coef = inp_flt(c_str(Diffiso_coef), diffiso_coefdef) diff_coef_tr = diffiso_coef diffiso_coef = inp_flt(c_str(Diffiso_coef_mo), diff_coef_tr) diff_coef_mo = diffiso_coef diff_dz_mo = inp_flt(c_str(Diffiso_dz_mo), 10000.) diffiso_slmax = inp_flt(c_str(Diffiso_slmax), diffiso_slmaxdef) slred=inp_flt(c_str(Diffiso_slred),diffiso_slreddef) sigzmin=inp_flt(c_str(Diffiso_sigzmin),sigzmindef) if (sigzmin.ge.0) then print*,'Diffiso_sigzmin should be negative' stop endif delta_rho = inp_flt(c_str(Diffiso_delta_rho), delta_rhodef) f0 = 0. tscl = inp_days (c_str(Diffiso_time), diffiso_tscldef) if (tscl.gt.delt/4.) then tscl = (delt/tscl)/4. else tscl = 1. endif diffiso_cadv=inp_flt(c_str(Diffiso_cadv),diff_coef_tr) c use_diff_cadv = diffiso_alpha.ne.0 .and. diffiso_cadv.ne.0 use_diff_cadv = diffiso_alpha.ne.0 if (use_diff_cadv) then facz_cnst = inp_flt(c_str(Diffiso_facz_cnst), facz_cnstdef) do k = 1, nz facz(k) = facz_cnst enddo nret = inp_rarr(c_str(Diffiso_facz),nz,facz,facz) psi_rel=inp_flt(c_str(Diffiso_psi_relax),psi_relaxdef) print*,'using Gent-McWilliams mixing ...' endif if (use_trdiff.ne.0) print*,'diff_coef_tr = ',diff_coef_tr if (use_modiff.ne.0) print*,'diff_coef_mo = ',diff_coef_mo endif diff_coef_dd = inp_flt(c_str(Div_damp_coef),0.) isteer = inp_flt(c_str(Div_steering),0) use_div_damp = diff_coef_dd.ne.0. use_tvd_vert = inp_def(c_str(TVD_vert)) use_tvd_horz = inp_def(c_str(TVD_horz)) cnst_upwind_z = inp_flt(c_str(Upwind_cnst_z),cnst_upwinddef) cnst_upwind_z_ts = inp_flt(c_str(Upwind_cnst_ts),cnst_upwind_z) cnst_upwind_z_tr = inp_flt(c_str(Upwind_cnst_tr),cnst_upwind_z) cnst_upwind_xy = inp_flt(c_str(Upwind_cnst_xy),cnst_upwinddef) cnst_upwind_xy_ts = inp_flt(c_str(Upwind_cnst_xy_ts),cnst_upwind_xy) cnst_upwind_xy_tr = inp_flt(c_str(Upwind_cnst_xy_tr),cnst_upwind_xy) if (cnst_upwind_z.lt.1. .or. cnst_upwind_xy.lt.1.or. * cnst_upwind_z_tr.lt.1. .or. cnst_upwind_xy_tr.lt.1.or. * cnst_upwind_z_ts.lt.1. .or. cnst_upwind_xy_ts.lt.1) then print*,'do not use Upwind_cnst less than one' stop endif cupi = 1./(cnst_upwind_z + 1.) cupi_ts = 1./(cnst_upwind_z_ts + 1.) cupi_tr = 1./(cnst_upwind_z_tr + 1.) hupi = 1./(cnst_upwind_xy + 1.) hupi_ts = 1./(cnst_upwind_xy_ts + 1.) hupi_tr = 1./(cnst_upwind_xy_tr + 1.) initb = inp_int(c_str(Bathymetry), 0) if (initb.ge.3) then i_ridge_min = inp_int(c_str(Ridge_min),0) i_ridge_max = inp_int(c_str(Ridge_max),nxp+1) endif initg = inp_int(c_str(Bath_grad), 2) grad_fac = inp_flt(c_str(Bath_grad_fac), 1.) initbt = inp_int(c_str(Bath_type), 0) depth_fac = inp_flt(c_str(Bath_fac), .2) nzz = inp_rarr(c_str(Z_profile), nz, zin, zin) nzh = inp_rarr(c_str(H_profile), nz, hin, hin) if (nz .gt. nzz .and. nz .gt. nzh) 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 do j = 1, nz if (hin(j).le.0) then print*,'a layer depth is negative, check Z_profile or H_profile' stop endif enddo 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 if (irest .eq. 0) then write(iout,*)'Z_Profile=',(zin(j),j=1,nz) write(iout,*)'H_Profile=',(hin(j),j=1,nz) write(iout,*)'W_Profile=',(zin(j)+dzin(j+1),j=1,nz) endif call flush (iout) dep = 0. do j = 1, nz dep = dep + hin(j) if (dep.le.tr_zmin) then ttr_fac(j)= 1. elseif (dep.ge.tr_zmax) then ttr_fac(j)= tr_fac else ttr_fac(j)= 1.+ (tr_fac - 1.)* (dep- tr_zmin) /(tr_zmax-tr_zmin) endif if (ttr_fac(j).ne.1.) print*,j,dep,ttr_fac(j) enddo trans_coef_sst = 0. trans_coef_sss = 0. if (rlx_time_sst.ne.0) trans_coef_sst = (zin(1)+zin(2))/2./rlx_time_sst if (rlx_time_sss.ne.0) trans_coef_sss = (zin(1)+zin(2))/2./rlx_time_sss print*,'Haney surface coefs, (T,S)',trans_coef_sst,trans_coef_sss t_in(nz+1) = inp_flt(c_str(Temp_bot), 0.) burger = inp_flt(c_str(Temp_burger), 0.) if (burger.ne.0.) then del_rho = - burger**2 * 0.14158 print*,'temperature profile:' do k = 1, nz z = zin(k) t_in(k) = t_in(nz+1) + 3.81847*(z-4500.)/4500. * del_rho print*,k,zin(k),t_in(k) enddo else i = inp_rarr(c_str(T_profile), nz, t_in, t_in) if (initt.le.1.and.i .lt. nz) then call perror1('T_profile: not enough terms...',-1) endif endif TATM = inp_flt(c_str(Temp_atm), 30.) i = inp_rarr(c_str(S_profile), nz, s_in, s_in) if (inits.le.1.and.i .lt. nz) then call perror1('S_profile: not enough terms...',-1) endif s_in(nz+1) = inp_flt(c_str(Salt_bot), 36.) 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(Pdens_bot),pdens_pnt(TEMP_BOT,SALT_BOT)) SITUD_BOT=inp_flt(c_str(Sdens_bot),sdens_pnt(TEMP_BOT,SALT_BOT,zin(nz+1))) rk_bot = inp_flt (c_str(Vert_kbot), rk_botdef) 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 nret=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') call mem_alloc (p_inx, nxyc, 1, 'inx') call mem_alloc (p_iny, nxyc, 1, 'iny') npt = nxyc npt1 = 1 npt2 = 1 + npt npt3 = 1 + 2*npt npt4 = 1 + 3*npt first_step = .true. call hfx_pert_init call wind_stir_init call nao_init ntrac_sur = 0 if (use_ice) then ntrac_ice = inp_int(c_str(Ice_ntrac),3) call ice_init(nx,ny,ntrac_ice,npt) ntrac_sur = ntrac_sur + ntrac_ice endif if (lev_err .ge. 1) then write (iout, *) 'NX =', nxp write (iout, *) 'NY =', nyp write (iout, *) 'NZ =', nz write (iout, *) 'NPACK =', nxyc write (iout, *) 'NTRACERS_KEITH =', ntrac write (iout, *) 'NTRACERS_BIO =', ntrac_bio 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(1000) character*465 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' * .or.ch.eq.'b'.or.ch.eq.'t') 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.'t') ich = 4 if (ch.eq.'b') ich = 5 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 c------------------------------------------------------------ subroutine bio_init c------------------------------------------------------------ include 'comm_para.h' include 'comm_new.h' include 'comm_data.h' include 'comm_bio.h' include 'biology.h' real inp_flt, inp_days logical inp_def dimension flt(100) c n_bio = inp_str(c_str(Nitrate_file) , c_str(nitrate_data), fbnitrate) JPBIO = inp_int(c_str(Bio_JPBIO),JPBIOdef) NFLBIO = inp_int(c_str(Bio_NFLBIO),NFLBIOdef) NDET = inp_int(c_str(Bio_NDET),NDETdef) NZOO = inp_int(c_str(Bio_NZOO),NZOOdef) NPHY = inp_int(c_str(Bio_NPHY),NPHYdef) NNUT = inp_int(c_str(Bio_NNUT),NNUTdef) APMIN = inp_flt(c_str(Bio_APMIN),APMINdef) AZMIN = inp_flt(c_str(Bio_AZMIN),AZMINdef) ANMIN = inp_flt(c_str(Bio_ANMIN),ANMINdef) ADMIN = inp_flt(c_str(Bio_ADMIN),ADMINdef) REDF = inp_flt(c_str(Bio_REDF),REDFdef) SLOPET = inp_flt(c_str(Bio_SLOPET),SLOPETdef) TOPTP = inp_flt(c_str(Bio_TOPTP),TOPTPdef) AKNUT = inp_flt(c_str(Bio_AKNUT),AKNUTdef) RCCHL = inp_flt(c_str(Bio_RCCHL),RCCHLdef) RGAMMA = inp_flt(c_str(Bio_RGAMMA),RGAMMAdef) TOPTGZ = inp_flt(c_str(Bio_TOPTGZ),TOPTGZdef) TMAXGZ = inp_flt(c_str(Bio_TMAXGZ),TMAXGZdef) RGZ = inp_flt(c_str(Bio_RGZ),RGZdef) RPPZ = inp_flt(c_str(Bio_RPPZ),RPPZdef) TAUS = inp_flt(c_str(Bio_TAUS),TAUSdef) AKS = inp_flt(c_str(Bio_AKS),AKSdef) FILMAX = inp_flt(c_str(Bio_FILMAX),FILMAXdef) RPNAZ = inp_flt(c_str(Bio_RPNAZ),RPNAZdef) RDNAZ = inp_flt(c_str(Bio_RDNAZ),RDNAZdef) EGGZOO = inp_flt(c_str(Bio_EGGZOO),EGGZOOdef) TAUZN = inp_flt(c_str(Bio_TAUZN),TAUZNdef) TMMAXP = inp_flt(c_str(Bio_TMMAXP),TMMAXPdef) TMMINP = inp_flt(c_str(Bio_TMMINP),TMMINPdef) TMMAXZ = inp_flt(c_str(Bio_TMMAXZ),TMMAXZdef) TMMINZ = inp_flt(c_str(Bio_TMMINZ),TMMINZdef) ANUMIN = inp_flt(c_str(Bio_ANUMIN),ANUMINdef) AFDMIN = inp_flt(c_str(Bio_AFDMIN),AFDMINdef) TAUDN = inp_flt(c_str(Bio_TAUDN),TAUDNdef) VSED = inp_flt(c_str(Bio_VSED),VSEDdef) AKI = inp_flt(c_str(Bio_AKI),AKIdef) TMAXR = inp_flt(c_str(Bio_TMAXR),TMAXRdef) TMINR = inp_flt(c_str(Bio_TMINR),TMINRdef) ALPH = inp_flt(c_str(Bio_ALPH),ALPHdef) TMUMAX = inp_flt(c_str(Bio_TMUMAX),TMUMAXdef) c c NAMOPT OPTICAL PARAMETERS c c XKG0 green absorption coefficient of water c XKR0 red absorption coefficent of water c XKGP green absorption coefficient of chl c XKRP red absorption coefficient of chl c XLG green chl exposant for absorption c XLR red chl exposant for absorption c RPIG chlahla+pheo ratio c XKD_A coefficient to take diffused light into account c XPAR2PUR_R coeff to transform par to pur in red c XPAR2PUR_G coeff to transform par to pur in green c XKG0 = inp_flt(c_str(Bio_XKG0),XKG0def) XKR0 = inp_flt(c_str(Bio_XKR0),XKR0def) XKGP = inp_flt(c_str(Bio_XKGP),XKGPdef) XKRP = inp_flt(c_str(Bio_XKRP),XKRPdef) XLG = inp_flt(c_str(Bio_XLG),XLGdef) XLR = inp_flt(c_str(Bio_XLR),XLRdef) RPIG = inp_flt(c_str(Bio_RPIG),RPIGdef) XKD_A = inp_flt(c_str(Bio_XKD_A),XKD_Adef) XPAR2PUR_R = inp_flt(c_str(Bio_XPAR2PUR_R),XPAR2PUR_Rdef) XPAR2PUR_G = inp_flt(c_str(Bio_XPAR2PUR_G),XPAR2PUR_Gdef) c c c RAU0 density of reference c RAU0 = inp_flt(c_str(Bio_RAU0),RAU0def) DELTAR = inp_flt(c_str(Bio_DELTAR),DELTARdef) return end c------------------------------------------------------------ subroutine ice_init(nx,ny,ntrac_ice,npt) c------------------------------------------------------------ include 'comm_new.h' include 'comm_data.h' include 'comm_para.h' include 'comm_amlice.h' include 'amlice.h' include 'dyice.h' include 'icedyn.h' real inp_flt, inp_days integer inp_rarr logical inp_def dimension flt(100) do i = 1, ntrac_ice ftrsnm_def(i) = "ICE_TR" enddo n_cice = 1 n_hice = 2 n_thice = 3 ftrsnm_def(n_cice) = "CICE" ftrsnm_def(n_hice) = "HICE" ftrsnm_def(n_thice)= "THICE" i_cice = npt*(n_cice - 1) + 1 i_hice = npt*(n_hice - 1) + 1 i_thice = npt*(n_thice- 1) + 1 iret= inp_sarr(c_str(Trs_name),ntrac_ice,ftrsnm_def,80,name_trs,ftrsnm) c default values are given amlice.h albedoocean = inp_flt(c_str(ICE_albedoocean), albedooceandef) albedoice = inp_flt(c_str(ICE_albedoice ), albedoicedef) albedof = inp_flt(c_str(ICE_albedof ), albedofdef) tfreeze = inp_flt(c_str(ICE_tfreeze ), tfreezedef) cicemax = inp_flt(c_str(ICE_cicemax ), cicemaxdef) hsnow = inp_flt(c_str(ICE_hsnow ), hsnowdef) sice = inp_flt(c_str(ICE_sice ), sicedef) itermax = inp_int(c_str(ICE_itermax ), itermaxdef) ssticemax = inp_flt(c_str(ICE_ssticemax ), ssticemaxdef) hicemin = inp_flt(c_str(ICE_hicemin ), hicemindef) tksnow = inp_flt(c_str(ICE_tksnow ), tksnowdef) tkice = inp_flt(c_str(ICE_tkice ), tkicedef) tkocean = inp_flt(c_str(ICE_tkocean ), tkoceandef) hq = inp_flt(c_str(ICE_hq ), hqdef) hf = inp_flt(c_str(ICE_hf ), hfdef) if (use_dyice) then dyice_lat = inp_flt(c_str(Dyice_lat),dyice_lat_def) Cdair = inp_flt(c_str(Dyice_Cdair),Cdair_def) Cdwater = inp_flt(c_str(Dyice_Cdwater),Cdwater_def) Pstar = inp_flt(c_str(Dyice_Pstar),Pstar_def) C = inp_flt(c_str(Dyice_C),C_def) phi = inp_flt(c_str(Dyice_phi),phi_def) delta = inp_flt(c_str(Dyice_delta),delta_def) etamax = inp_flt(c_str(Dyice_etamax),etamax_def) theta_a = inp_flt(c_str(Dyice_theta_a),theta_a_def) theta_w = inp_flt(c_str(Dyice_theta_w),theta_w_def) do k = 1, ntrac_ice Ktracer(k) = 0.0 ! diff coeff for h (ice thick) [m2/s] enddo nic = inp_rarr(c_str(Dyice_Ktracer),0,Ktracer,Ktracer) BndyCond = inp_int(c_str(Dyice_BndyCond),1) DragLaw = inp_int(c_str(Dyice_DragLaw),2) Rheology = inp_int(c_str(Dyice_Rheology),3) InitialGuess = inp_int(c_str(Dyice_InitialGuess),2) endif if (use_dyice_old) then dyice_diff = inp_flt(c_str(Ice_dyice_diff),dyice_diffdef) endif return end c--------------------------------------------------------- subroutine baro_init c--------------------------------------------------------- include 'comm_para.h' include 'comm_new.h' include 'comm_data.h' 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 /errors/ ioerr, nstep 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_DEPMIN, * BAR_DSINK, ibar_key, nbaro, rayl, nonlin_baro character*1 BC_I(0:9) common /masks/BC_I logical use_per_island, use_stan_island, use_sunk common /baro_island/ * alons_min(10),alons_max(10),alats_min(10),alats_max(10) * ,alon1_min,alon1_max,alat1_min,alat1_max * ,per_lat,use_sunk,use_per_island,use_stan_island * ,itrans,ftrans common /new_island/n_islands,in_rhs(0:9),ityp_island(0:9), * alon_island(9),alat_island(9),btrans(0:9), * n_poly(9),alon_poly(100,9),alat_poly(100,9), * i_poly(100,9),j_poly(100,9) c-------------------------------------------------------------- logical sink_region common /bathymetry/sink_region,sr_depth, alon_sr, blon_sr, alat_sr, blat_sr common/friction/ ibfric,b_fric, sw_fric integer inp_rarr real inp_flt, inp_days, xf(100) logical inp_def character*10 numbers numbers = '1234567890' 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) nonlin_baro= inp_int(c_str(Baro_nonlin), 1) if (nonlin_baro.ne.1) then print*,'Baro_nonlin 0 no longer allowed' stop endif BAR_DSINK = inp_flt(c_str(Baro_depsink), dep_min - 1.) BAR_DEPMIN = inp_flt(c_str(Baro_depmin), dep_min) rayl = inp_flt(c_str(Baro_rayl), 0.) ibar_key = inp_int(c_str(Baro_err), 1) adelt = inp_flt(c_str(Island_relax), 0.) if (ibar_key .ne. 0) * n_bar = inp_str(c_str(Baro_errfile), * finp(1:mlen(finp))//'.bar\0', f_bar) use_per_island = inp_def(c_str(Antarctica)) .or. * inp_def(c_str(Channel)) if (use_per_island) then ityp_island(0) = 1 ! no sinking of periodic island allowed if (iglob.eq.0) then print*,'cannot have periodic island in non-periodic domain,' print*,' - check Periodic and Periodic_island settings' stop endif if (inp_def(c_str(Antarctica))) then per_lat = inp_flt(c_str(Antarctica_lat),-60) atrans = inp_flt(c_str(Antarctica_transport),0.) if (atrans.ne.0) ityp_island(0) = 3 elseif (inp_def(c_str(Channel))) then per_lat = inp_flt(c_str(Channel_lat),(alat+blat)/2.) atrans = inp_flt(c_str(Channel_transport),0.) itrans = inp_int(c_str(Channel_transport_type),0) ftrans = inp_days(c_str(Channel_transport_freq),1.) if (atrans.ne.0) ityp_island(0) = 3 endif btrans(0) = atrans endif nis = inp_rarr(c_str(New_island),0,aflag,aflag) n_islands = aflag(1) if (n_islands.gt.9) then print*,'since you are using more than 9 islands' print*,'you must allocate more islands in barotropic.h' stop endif if (n_islands.gt.0.and.nis.ne.3*n_islands+1) then print*,'you must specify a long,lat and island type for each island' stop endif use_sunk = .false. use_stan_island = .false. do kis=1,n_islands ityp_island(kis) = aflag(3*(kis-1)+2) if (ityp_island(kis).eq.2) use_sunk = .true. if (ityp_island(kis).eq.1.or.ityp_island(kis).eq.3) * use_stan_island = .true. alon_island(kis) = aflag(3*(kis-1)+3) alat_island(kis) = aflag(3*(kis-1)+4) BC_I(kis) = numbers(kis:kis) nis = inp_rarr(n_a(Island_poly,kis),0,xf,xf) btrans(kis) = inp_flt(n_a(Island_tran,kis),0.) n_poly(kis) = nis/2 if (nis/2 .ne. float(nis)/2.) then print*,'check long,lat points for New_island_poly',kis stop endif do i = 1, n_poly(kis) alon_poly(i,kis) = xf(2*(i-1)+1) alat_poly(i,kis) = xf(2*i) enddo enddo sink_region = .false. nsr = inp_rarr(c_str(Sink_region),0,aflag,aflag) if (nsr.gt.0) then sink_region = .true. if (nsr.ne.5) then print*,'For Sink_region, must specify * long(min,max),lat(min,max),depth' stop endif sr_depth = aflag(5) alon_sr = aflag(1) blon_sr = aflag(2) alat_sr = aflag(3) blat_sr = aflag(4) endif 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. nret = inp_iarr (c_str(Baro_iflag), 5, iflag, iflag) nret = inp_rarr (c_str(Baro_aflag), 5, aflag, aflag) return end