c ------------------------------------------------ subroutine depth_init (npt, zin) c ------------------------------------------------ implicit real(a-h,o-z),integer(i-n) include 'comm_data.h' include 'comm_new.h' 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 dimension tp(1), zin(1) if (initb .eq. 0) then c.....Constant Depth (Flat Bottom) dep_max = zin(nz+1) do i = 1, npt dept(i) = dep_max enddo dep_min = dep_max elseif (initb .eq. 1) then c.....Read Bathymetry Data From a File call odb_open(idf_dp, fbdep(1:n_dep), 0) call data_on_model_grid(idf_dp, lret, 'bath') call read_zt (idf_dp, lret, npt, 1, 1, 'bath', tp, dept) c do i = 1, npt c dept(i) = max(dept(i),(zin(3)+zin(2))/2.) c dept(i) = max(dept(i),(zin(4)+zin(5))/2.) c enddo elseif (initb .eq. 2) then dept(1) = 50 + zin(nz+1)*ran(10001) do i = 2, npt dept(i) = 50 + zin(nz+1)*ran( ) enddo do i = 1, npt dept(i) = max(dept(i),zin(1)) enddo elseif (initb .ge. 3) then c ramp depmin = (zin(initb-2) + zin(initb-1))/2. depmax = zin(nz+1) x1 = xm(1) xn = xm(nxp) xs = xn-x1 do i = 1, nxp tp(i) = depmin + (xm(i) - x1)*(depmax-depmin)/xs enddo do k = 1, npt j = (iox(k)-1)/nxp + 1 i = iox(k) - (j-1)*nxp dept(k) = tp(i) enddo endif dept(1) = max(dep_min,dept(1)) do i = 2, npt dept(i) = max(dep_min,dept(i)) enddo dep_min = dept(1) dep_max = dept(1) do i = 2, npt dep_max = max(dep_max, dept(i)) dep_min = min(dep_min, dept(i)) enddo return end c -------------------------------------------------------------------- subroutine clim_init(npt,nstart,h0,sigma,dzin,hmf, * hclf,tclf,sclf,dclf,tpf,nsponge,lsponge) c -------------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' include 'comm_data.h' dimension h0(1),hmf(npt,nz),hclf(npt,nz,1),dclf(npt,nz),dzin(nz+1), * tclf(npt,nz,1),sclf(npt,nz,1),tpf(npt,1),ind3(3),sigma(nz) * ,lsponge(nsponge) 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 c do k = nsig+1, nz do i = 1, npt hclf(i,k,1) = hmf(i,k) enddo enddo if (icl_h .eq. 0) then do k = 1, nsig do i = 1, npt hclf(i,k,1) = hmf(i,k) enddo enddo elseif (icl_h .eq. 1) then c........H: sigma structure between Mixed Layer & Thermocline - Static: call odb_open(idf_hcl, fbhcl(1:n_hcl), 0) call data_on_model_grid(idf_hcl, lclm, 'mltc') if (icl_htop .eq. 1) * call read_zt (idf_hcl,lclm,npt, 1,1, 'mltc', tpf, hclf(1,1,1)) if (icl_htop .eq. 0) call afill (npt, hclf, h0(1)) sigk = sigma(3) do i = 1, npt hclf(i,2,1) = 0.5*hclf(i,1,1) + sigk*(z_begin - 1.5*hclf(i,1,1)) enddo do k = 3, nsig - 1 sigkp = sigma(k+1) do i = 1, npt hclf(i,k,1) = (sigk+sigkp) * (z_begin - 1.5*hclf(i,1,1)) enddo sigk = sigkp enddo k = nsig do i = 1, npt hclf(i,k,1) = dzin(k+1) + sigma(k)*(z_begin - 1.5*hclf(i,1,1)) enddo elseif (icl_h .eq. 2) then c........H: sigma structure between Mixed Layer & Thermocline - Dynamic call odb_open(idf_hcl, fbhcl(1:n_hcl), 0) call data_on_model_grid(idf_hcl, lclm, 'mltc') call odb_rddm(idf_hcl, 'T', ntclm) call mem_alloc(p_tclm, ntclm, 2, 'tclm') call odb_rdgr(idf_hcl, 'T', ntclm, tclm) call it_catch (ntclm, tclm, nstart, it1, it2, clm_tscl) iclm = it2 if (icl_htop .eq. 1) then call read_zt (idf_hcl, lclm,npt, 1, it1, 'mltc', tpf, hclf(1,1,1)) call read_zt (idf_hcl, lclm,npt, 1, it2, 'mltc', tpf, hclf(1,1,2)) endif if (icl_htop .eq. 0) then call afill (npt, hclf(1,1,1), h0(1)) call afill (npt, hclf(1,1,2), h0(1)) endif sigk = sigma(3) do i = 1, npt hclf(i,2,1) = 0.5*hclf(i,1,1) + sigk*(z_begin - 1.5*hclf(i,1,1)) hclf(i,2,2) = 0.5*hclf(i,1,2) + sigk*(z_begin - 1.5*hclf(i,1,2)) enddo do k = 3, nsig - 1 sigkp = sigma(k+1) do i = 1, npt hclf(i,k,1) = (sigk+sigkp) * (z_begin - 1.5*hclf(i,1,1)) hclf(i,k,2) = (sigk+sigkp) * (z_begin - 1.5*hclf(i,1,2)) enddo sigk = sigkp enddo k = nsig do i = 1, npt hclf(i,k,1) = dzin(k+1) + sigma(k)*(z_begin - 1.5*hclf(i,1,1)) hclf(i,k,2) = dzin(k+1) + sigma(k)*(z_begin - 1.5*hclf(i,1,2)) enddo endif if (icl_ts .eq. 0) return call odb_open(idf_t, fbtem(1:n_tem), 0) if (icl_h .eq. 0) call data_on_model_grid(idf_t, lclm, 'temp') call odb_rddm(idf_t, 'Z', nzclm) call mem_alloc(p_zclm, nzclm, 2, 'MZ for temp climatology') call odb_rdgr (idf_t, 'Z', nzclm, zclm) mz = nzclm if (use_salt) then call odb_open(idf_s, fbsal(1:n_sal), 0) call odb_rddm(idf_s, 'Z', mz) if (mz .ne. nzclm) * call perror1('Temp & Salt DATA should be on the same Z grid', 1) endif if (icl_ts .eq. 1) then !!.....time independent case: call read_linz(idf_t,lclm,npt,mpack,nz,mz,1,hclf,tclf,zclm,tp,'temp') if (use_salt) * call read_linz(idf_s,lclm,npt,mpack,nz,mz,1,hclf,sclf,zclm,tp,'salt') if (iv_bot .eq. 4) then do i = 1, npt tclf(i,nz,1) = TEMP_BOT if (use_salt) sclf(i,nz,1) = SALT_BOT enddo endif elseif (icl_ts .eq. 2) then !!.....time dependent case: call odb_rddm(idf_t, 'T', i) if (icl_h .ne. 2) then !!.....H_clim is time independent ntclm = i call mem_alloc(p_tclm, ntclm, 2, 'tclm') call odb_rdgr(idf_t, 'T', ntclm, tclm) call it_catch (ntclm, tclm, nstart, it1, it2, clm_tscl) iclm = it2 elseif (i .ne. ntclm ) then call perror1('MxTh & Temp DATA should be on the same T grid', 1) endif if (use_salt) then call odb_rddm(idf_s, 'T', i) if (i .ne. ntclm) * call perror1('Temp & Salt DATA should be on the same T grid', 1) endif call read_linz(idf_t, lclm, npt,mpack,nz,mz, it1, * hclf(1,1,1),tclf(1,1,1),zclm,tp, 'temp') k2_h = 1 if (icl_h .eq. 2) k2_h = 2 call read_linz(idf_t, lclm, npt,mpack,nz,mz, it2, * hclf(1,1,k2_h),tclf(1,1,2),zclm,tp, 'temp') if (use_salt) then call read_linz(idf_s, lclm, npt,mpack,nz,mz, it1, * hclf(1,1,1),sclf(1,1,1),zclm,tp, 'salt') call read_linz(idf_s, lclm, npt,mpack,nz,mz, it2, * hclf(1,1,k2_h),sclf(1,1,2),zclm,tp, 'salt') endif endif if (ipre.eq.1) then if (use_salt) call potn_dens (npt,nzi,tclf(1,1,1),sclf(1,1,1),dclf) call dconv_cl (npt,nz,nzi,hclf,tclf(1,1,1),sclf(1,1,1),dclf) if (use_salt) call potn_dens (npt,nzi,tclf(1,1,2),sclf(1,1,2),dclf) call dconv_cl (npt,nz,nzi,hclf,tclf(1,1,2),sclf(1,1,2),dclf) endif if (icl_rlx .eq. 1) then do j = 1, nyp if (ym(j) .gt. clm_no) then tp(j) = clm_coef * (ym(j)-clm_no)/(ym(nyp)-clm_no) elseif (ym(j) .lt. clm_so) then tp(j) = clm_coef * (clm_so-ym(j))/(clm_so-ym(1)) else tp(j) = 0. endif enddo do i = 1, npt j = (iox(i)-1)/nxp + 1 relax(i) = tp(j) enddo elseif (icl_rlx .eq. 2) then do i = 1, npt relax(i) = clm_coef enddo elseif (icl_rlx .eq. 3) then dsponge = real(ksponge) do i = 1, npt ixy = iox(i) ix = mod (ixy -1 ,nxp) + 1 iy = (ixy - ix)/nxp + 1 dmin = nxp+nyp do j = 1, nsponge kxy = lsponge(j) kx = mod (kxy -1 ,nxp) + 1 ky = (kxy - kx)/nxp + 1 d = abs(ix-kx) + abs(iy-ky) dmin = min(d,dmin) enddo relax(i) = clm_coef*max((dsponge-dmin)/dsponge,0.) enddo endif return end c ---------------------------------------------- subroutine h_init (npt, nz, nstart, hmf, hclf) c ---------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' dimension hmf(npt,1), hclf(npt,nz,1) if (icl_h .eq. 1) then do k = 1, nz do i = 1, npt hmf(i,k) = hclf(i,k,1) enddo enddo elseif (icl_h .eq. 2) then do k = 1, nz do i = 1, npt hmf(i,k) = hclf(i,k,1) + clm_tscl*(hclf(i,k,2) - hclf(i,k,1)) enddo enddo endif return end c ----------------------------------------------------------------- subroutine temp_init (npt, nz, nzi, nstart, t0, tmf, tclf) c ----------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' dimension t0(1), tmf(npt,1), tclf(npt,nz,1), nzi(1) c if (initt .eq. 0) then c.......constant values for each layer according to T_INIT do i = 1, npt do k = 1, nzi(i) tmf(i,k) = t0(k) enddo enddo elseif (initt .eq. 3) then c.....from CLIMATOLOGY Data. if (icl_ts .eq. 2) then do i = 1, npt do k = 1, nzi(i) tmf(i,k) = tclf(i,k,1)+ clm_tscl*(tclf(i,k,2)- tclf(i,k,1)) enddo enddo else do i = 1, npt do k = 1, nzi(i) tmf(i,k) = tclf(i,k,1) enddo enddo endif endif return end c ----------------------------------------------------------------- subroutine salt_init (npt, nz, nzi, nstart, s0, smf, sclf) c ----------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' dimension s0(1), smf(npt,1), sclf(npt,nz,1), nzi(1) c if (inits .eq. 0) then do i = 1, npt do k = 1, nzi(i) smf(i,k) = s0(k) enddo enddo elseif (inits .eq. 3) then c.....from CLIMATOLOGY Data. if (icl_ts .eq. 2) then do i = 1, npt do k = 1, nzi(i) smf(i,k) = sclf(i,k,1)+ clm_tscl*(sclf(i,k,2)- sclf(i,k,1)) enddo enddo else do i = 1, npt do k = 1, nzi(i) smf(i,k) = sclf(i,k,1) enddo enddo endif endif return end c---------------------------------------------------- subroutine tau_init (nstart,npt, dtaux, dtauy) c---------------------------------------------------- c..........initalize the winds according to MTAU implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' include 'comm_data.h' parameter (TAUCON = 10300., API = 3.14159265, TAUINV = 1./TAUCON) common/winds/mtau,matau,tausc,atau,froude,icloud 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 dimension dtaux(npt,1),dtauy(npt,1) taub = tausc/TAUCON tauc = atau/TAUCON if (mtau .eq. 1 .or. mtau .eq. 2) then call odb_open(idf_tx, fbwnd(1:n_wnd)//'.x', 0) call odb_open(idf_ty, fbwnd(1:n_wnd)//'.y', 0) call odb_rddm(idf_tx, 'T', ntau) call mem_alloc(p_ttau, ntau, 2, 'ttau') call odb_rdgr(idf_tx, 'T', ntau, ttau) call it_catch (ntau, ttau, nstart, it1, it2, tscl) itau = it2 call data_on_model_grid(idf_tx, ltau, 'tau') call read_zt (idf_tx, ltau, npt, 1, it1, 'taux', tp, dtaux(1,1)) call read_zt (idf_tx, ltau, npt, 1, it2, 'taux', tp, dtaux(1,2)) call read_zt (idf_ty, ltau, npt, 1, it1, 'tauy', tp, dtauy(1,1)) call read_zt (idf_ty, ltau, npt, 1, it2, 'tauy', tp, dtauy(1,2)) do i = 1, npt dtaux(i,1) = TAUINV * dtaux(i,1) dtaux(i,2) = TAUINV * dtaux(i,2) dtauy(i,1) = TAUINV * dtauy(i,1) dtauy(i,2) = TAUINV * dtauy(i,2) enddo if (mtau .eq. 1) then do i = 1, npt taux(i) = dtaux(i,1) + tscl * (dtaux(i,2) - dtaux(i,1)) tauy(i) = dtauy(i,1) + tscl * (dtauy(i,2) - dtauy(i,1)) enddo endif elseif (mtau .eq. 3) then c..........3 - annualy averaged climatology call odb_open(idf_tx, fbwnd(1:n_wnd)//'.x', 0) call odb_open(idf_ty, fbwnd(1:n_wnd)//'.y', 0) call odb_rddm(idf_tx, 'T', ntau) call data_on_model_grid(idf_tx, ltau, 'tau') do k = 1, ntau call read_zt (idf_tx, ltau, npt, 1, k, 'taux', tp, dtaux(1,2)) call read_zt (idf_ty, ltau, npt, 1, k, 'tauy', tp, dtauy(1,2)) do i = 1, npt dtaux(i,1) = dtaux(i,1) + dtaux(i,2) dtauy(i,1) = dtauy(i,1) + dtauy(i,2) enddo enddo coef = TAUINV/real(ntau) do i = 1, npt taux(i) = coef * dtaux(i,1) tauy(i) = coef * dtauy(i,1) enddo elseif (mtau .eq. 5) then c..........5 - COSINE winds if (itau_cos .eq. 0) then do j = 1, nyp tp(j) = taub*cos(API*(ym(j))/80.) enddo elseif (itau_cos .eq. 1) then y1 = ym(1) y2 = ym(nyp) do j = 1, nyp tp(j) = taub*cos(2.*API*( (ym(j)-y1)/(y2-y1) - 0.5)) enddo endif do k = 1, npt j = (iox(k)-1)/nxp + 1 tmpx = tp(j) taux(k) = tmpx dtaux(k,1) = tmpx tauy(k) = tauc dtauy(k,1) = tauc enddo endif return end c ------------------------------------------------------------------ subroutine tau_lin (nstep,npt,ixd,im2d,blcf, taux,tauy,dtx,dty,tp) c ------------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' parameter (TAUCON = 10300., API = 3.14159265, TAUINV = 1./TAUCON) dimension taux(1),tauy(1),dtx(npt,1),dty(npt,1) dimension ixd(1),im2d(1),blcf(1),tp(1) call it_catch (ntau, ttau, nstep, it1, it2, tscl) if (it2 .ne. itau) then itau = it2 do i = 1, npt dtx(i,1) = dtx(i,2) dty(i,1) = dty(i,2) enddo call read_zt (idf_tx, ltau, npt, 1, it2, 'taux', tp, dtx(1,2)) call read_zt (idf_ty, ltau, npt, 1, it2, 'tauy', tp, dty(1,2)) do i = 1, npt dtx(i,2) = TAUINV * dtx(i,2) dty(i,2) = TAUINV * dty(i,2) enddo endif do i = 1, npt taux(i) = dtx(i,1) + tscl * (dtx(i,2) - dtx(i,1)) tauy(i) = dty(i,1) + tscl * (dty(i,2) - dty(i,1)) enddo return end c ------------------------------------------------------------------ subroutine hflx_init (nstart, npt, nx, ny, temp, sstf, cldf, slrf) c ------------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' include 'comm_data.h' include 'comm_pbl.h' 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 dimension temp(1), sstf(npt,1), cldf(npt,1), slrf(npt,1) if (initq .eq. 0) then c.....Haney case Q = QCOF * (TATM - T) with Heat Transfer Coef.=30 w/m**2/k. do i = 1, npt sstf(i,1) = TATM enddo return elseif (initq .eq. 1) then c.....T_atm = initial T(i,1) = const(time) do i = 1, npt sstf(i,1) = temp(i) enddo elseif (initq .eq. 2) then c.....T_atm = average(SST) = const(time) call odb_open(idf_sst, fbsst(1:n_sst), 0) call odb_rddm(idf_sst, 'T', nsst) call data_on_model_grid(idf_sst, lsst, 'sst') do k = 1, nsst call read_zt (idf_sst, lsst, npt, 1, k, 'sst', tp, sstf(1,2)) do i = 1, npt sstf(i,1) = sstf(i,1) + sstf(i,2) enddo enddo coef = 1./real(nsst) do i = 1, npt sstf(i,1) = coef * sstf(i,1) enddo elseif (initq .eq. 3) then c.....T_atm = SST(time) - climatology call odb_open(idf_sst, fbsst(1:n_sst), 0) call odb_rddm(idf_sst, 'T', nsst) call mem_alloc(p_tsst, nsst, 2, 'tsst') call odb_rdgr(idf_sst, 'T', nsst, tsst) call it_catch (nsst, tsst, nstart, it1, it2, tscl) isst = it2 call data_on_model_grid(idf_sst, lsst, 'sst') call read_zt (idf_sst, lsst, npt, 1, it1, 'sst', tp, sstf(1,1)) call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', tp, sstf(1,2)) elseif (initq .eq. 4) then c.....Q = Q(CLD,TAU,SLR),using Seager formulation.+[qcof*(SST_clim - SST_model)] call odb_open(idf_sst, fbsst(1:n_sst), 0) call odb_open(idf_cld, fbcld(1:n_cld), 0) call odb_rddm(idf_sst, 'T', nsst) call mem_alloc(p_tsst, nsst, 2, 'tsst') call odb_rdgr(idf_sst, 'T', nsst, tsst) call odb_rddm(idf_cld, 'T', i) if (i .ne. nsst) * call perror1('SST & Cloud Cover DATA should be on the same grid',1) call it_catch (nsst, tsst, nstart, it1, it2, tscl) isst = it2 call data_on_model_grid(idf_sst, lsst, 'sst') call read_zt (idf_sst, lsst, npt, 1, it1, 'sst', tp, sstf(1,1)) call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', tp, sstf(1,2)) call read_zt (idf_cld, lsst, npt, 1, it1, 'cld', tp, cldf(1,1)) call read_zt (idf_cld, lsst, npt, 1, it2, 'cld', tp, cldf(1,2)) elseif (initq .eq. 5) then c.....*second* Richard-Benno formulation: Q = Q(T, solr, wndsp, clouds) call odb_open(idf_sst, fbsst(1:n_sst), 0) call odb_open(idf_cld, fbcld(1:n_cld), 0) call odb_open(idf_slr, fbslr(1:n_slr), 0) call odb_rddm(idf_sst, 'T', nsst) call mem_alloc(p_tsst, nsst, 2, 'tsst') call odb_rdgr(idf_sst, 'T', nsst, tsst) call odb_rddm(idf_cld, 'T', i) if (i .ne. nsst) * call perror1('H.flx & cloud data should be on the same grid',1) call odb_rddm(idf_slr, 'T', i) if (i .ne. nsst) * call perror1('H.flx & Sol.Rad. data should be on the same grid',1) call it_catch (nsst, tsst, nstart, it1, it2, tscl) isst = it2 call data_on_model_grid(idf_sst, lsst, 'sst') call read_zt (idf_sst, lsst, npt, 1, it1, 'sst', tp, sstf(1,1)) call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', tp, sstf(1,2)) call read_zt (idf_cld, lsst, npt, 1, it1, 'cld', tp, cldf(1,1)) call read_zt (idf_cld, lsst, npt, 1, it2, 'cld', tp, cldf(1,2)) call read_zt (idf_slr, lsst, npt, 1, it1, 'solr',tp, slrf(1,1)) call read_zt (idf_slr, lsst, npt, 1, it2, 'solr',tp, slrf(1,2)) elseif (initq .eq. 6) then trans_coef = 0. elseif (initq .eq. 7) then c.....SST = Annual Mean SST; Cld = Annual Mean Cloud Cover: call odb_open(idf_sst, fbsst(1:n_sst), 0) call odb_open(idf_cld, fbcld(1:n_cld), 0) call odb_rddm(idf_cld, 'T', nsst) call odb_rddm(idf_cld, 'T', ncld) if (nsst .ne. ncld) * call perror1('SST & Cloud Cover DATA should have same grids', 1) call data_on_model_grid(idf_sst, lsst, 'sst') do k = 1, nsst call read_zt (idf_sst, lsst, npt, 1, k, 'sst', tp, sstf(1,2)) call read_zt (idf_cld, lsst, npt, 1, k, 'cld', tp, cldf(1,2)) do i = 1, npt sstf(i,1) = sstf(i,1) + sstf(i,2) cldf(i,1) = cldf(i,1) + cldf(i,2) enddo enddo coef = 1./real(nsst) do i = 1, npt sstf(i,1) = coef * sstf(i,1) cldf(i,1) = coef * cldf(i,1) enddo elseif (initq .eq. 8) then c.....PBL model Q = Q(T, solr, wndsp, clouds) call odb_open(idf_sst, fbsst(1:n_sst), 0) call odb_open(idf_cld, fbcld(1:n_cld), 0) call odb_open(idf_slr, fbslr(1:n_slr), 0) call odb_open(idf_wsp, fwsp(1:n_wsp), 0) call odb_open(idf_uwd, fuwd(1:n_uwd), 0) call odb_open(idf_vwd, fvwd(1:n_vwd), 0) call odb_open(idf_ah, fah(1:n_ah), 0) call odb_open(idf_at, fat(1:n_at), 0) call odb_rddm(idf_sst, 'T', nsst) call odb_rddm(idf_cld, 'T', ncld) call odb_rddm(idf_slr, 'T', nslr) call odb_rddm(idf_wsp, 'T', nwsp) call odb_rddm(idf_uwd, 'T', nuwd) call odb_rddm(idf_vwd, 'T', nvwd) call odb_rddm(idf_ah, 'T', nah) call odb_rddm(idf_at, 'T', nat) if (nsst+ncld+nslr+nwsp+nuwd+nvwd+nah+nat .ne. 8*nsst) * call perror1('All PBL data should be on the same grid',1) call mem_alloc(p_tsst, nsst, 2, 'tsst') call odb_rdgr(idf_sst, 'T', nsst, tsst) call it_catch (nsst, tsst, nstart, it1, it2, tscl) isst = it2 call data_on_model_grid(idf_sst, lsst, 'sst') call read_zt (idf_sst, lsst, npt, 1, it1, 'sst', tp, sstf(1,1)) call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', tp, sstf(1,2)) call read_zt (idf_cld, lsst, npt, 1, it1, 'cld', tp, cldf(1,1)) call read_zt (idf_cld, lsst, npt, 1, it2, 'cld', tp, cldf(1,2)) call read_zt (idf_slr, lsst, npt, 1, it1, 'solr',tp, slrf(1,1)) call read_zt (idf_slr, lsst, npt, 1, it2, 'solr',tp, slrf(1,2)) call read_zt (idf_wsp, lsst, npt, 1, it1, 'wndspd', tp, wnsp(1,1)) call read_zt (idf_wsp, lsst, npt, 1, it2, 'wndspd', tp, wnsp(1,2)) call read_zt (idf_uwd, lsst, npt, 1, it1, 'uwnd', tp, uwnd(1,1)) call read_zt (idf_uwd, lsst, npt, 1, it2, 'uwnd', tp, uwnd(1,2)) call read_zt (idf_vwd, lsst, npt, 1, it1, 'vwnd', tp, vwnd(1,1)) call read_zt (idf_vwd, lsst, npt, 1, it2, 'vwnd', tp, vwnd(1,2)) NXY = nxp * nyp #ifdef gidon call read_zt (idf_ah, 0, NXY, 1, it1, 'spechum', tp, ahum(1,1)) call read_zt (idf_ah, 0, NXY, 1, it2, 'spechum', tp, ahum(1,2)) call read_zt (idf_at, 0, NXY, 1, it1, 'Tair', tp, atem(1,1)) call read_zt (idf_at, 0, NXY, 1, it2, 'Tair', tp, atem(1,2)) #else call read_zt (idf_ah, 0, NXY, 1, it1, 'airhum', tp, ahum(1,1)) call read_zt (idf_ah, 0, NXY, 1, it2, 'airhum', tp, ahum(1,2)) call read_zt (idf_at, 0, NXY, 1, it1, 'airtem', tp, atem(1,1)) call read_zt (idf_at, 0, NXY, 1, it2, 'airtem', tp, atem(1,2)) #endif endif return end c --------------------------------------------------- subroutine evpr_init (nstart, npt, salt, sssf, evpf) c --------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' include 'comm_data.h' 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 dimension salt(1), sssf(npt,1), evpf(npt,1) if (initep .eq. 0) then c given by EP = QCOF * (SATM - S) with Salt Transfer Coef.=30 o.e./m**2. do i = 1, npt sssf(i,1) = SATM enddo return elseif (initep .eq. 1) then c.....S_atm = initial S(i,1) = const(time) do i = 1, npt sssf(i,1) = salt(i) enddo elseif (initep .eq. 2) then c.....S_atm = average(SSS) = const(time) call odb_open(idf_sss, fbsss(1:n_sss), 0) call odb_rddm(idf_sss, 'T', nsss) if (nsst .ne. nsss) * call perror1('SST & SSS DATA should be on the same T grid',1) do k = 1, nsst call read_zt (idf_sss, lsst, npt, 1, k, 'sss', tp, sssf(1,2)) do i = 1, npt sssf(i,1) = sssf(i,1) + sssf(i,2) enddo enddo coef = 1./real(nsst) do i = 1, npt sssf(i,1) = coef * sssf(i,1) enddo elseif (initep .eq. 3) then c.....S_atm = SSS(time) - climatology call odb_open(idf_sss, fbsss(1:n_sss), 0) call odb_rddm(idf_sss, 'T', nsss) if (nsst .ne. nsss) * call perror1('SST & SSS DATA should be on the same T grid',1) call it_catch (nsst, tsst, nstart, it1, it2, tscl) isss = it2 call read_zt (idf_sss, lsst, npt, 1, it1, 'sss', tp, sssf(1,1)) call read_zt (idf_sss, lsst, npt, 1, it2, 'sss', tp, sssf(1,2)) elseif (initep .eq. 4) then c.....read EP directly: call odb_open(idf_evp, fbevp(1:n_evp), 0) call odb_rddm(idf_evp, 'T', nevp) call mem_alloc(p_tevp, nevp, 2, 'evp') call odb_rdgr(idf_evp, 'T', nevp, tevp) call it_catch (nevp, tevp, nstart, it1, it2, tscl) ievp = it2 call data_on_model_grid(idf_evp, levp, 'evp') call read_zt (idf_evp, levp, npt, 1, it1, 'evp', tp, evpf(1,1)) call read_zt (idf_evp, levp, npt, 1, it2, 'evp', tp, evpf(1,2)) elseif (initep .eq. 6) then trans_coef = 0. elseif (initep.eq.8 .or. initep.eq.9) then call odb_open(idf_evp, fbevp(1:n_evp), 0) call odb_rddm(idf_evp, 'T', nevp) if(nevp.ne.nsst) call perror1('evp should have PBL T grid',1) call it_catch (nsst, tsst, nstart, it1, it2, tscl) ievp = it2 call read_zt (idf_evp, lsst,npt,1,it1,'precip' ,tp, evpf(1,2)) call read_zt (idf_evp, lsst,npt,1,it2,'precip' ,tp, evpf(1,3)) endif return end c-------------------------------------------------------------------- subroutine qforc(nstep, npt, nx, ny, sstf, cldf, slrf, tpf, qbf) c--------------------------------------------------------------------- c update heat flux using current t(1), and SST(i) implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' include 'comm_data.h' include 'comm_pbl.h' dimension slrf(npt,1),sstf(npt,1),cldf(npt,1), tpf(npt,1), qbf(npt,1) if (initq .eq. 0) then do i = 1, npt q(i) = trans_coef * (TATM - t(i)) enddo elseif (initq .eq. 1 .or. initq .eq. 2) then do i = 1, npt q(i) = trans_coef * (sstf(i,1) - t(i)) enddo elseif (initq .eq. 3) then call it_catch (nsst, tsst, nstep, it1, it2, sst_tscl) if (it2 .ne. isst) then isst = it2 do i = 1, npt sstf(i,1) = sstf(i,2) enddo call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', tp, sstf(1,2)) endif do i = 1, npt sst_val = sstf(i,1) + sst_tscl * (sstf(i,2) - sstf(i,1)) sstf(i,3) = sst_val q(i) = trans_coef * (sst_val - t(i)) enddo elseif (initq .eq. 4) then call it_catch (nsst, tsst, nstep, it1, it2, tscl) if (it2 .ne. isst) then isst = it2 do i = 1, npt sstf(i,1) = sstf(i,2) cldf(i,1) = cldf(i,2) enddo call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', tp, sstf(1,2)) call read_zt (idf_cld, lsst, npt, 1, it2, 'cld', tp, cldf(1,2)) endif do i = 1, npt tpf(i,1) = sstf(i,1) + tscl * (sstf(i,2) - sstf(i,1)) tpf(i,2) = cldf(i,1) + tscl * (cldf(i,2) - cldf(i,1)) enddo tenso = enso_start + enso_scale * nstep call hflx_s89(tenso,npt,iox,t,tpf,tpf(1,2),ym, * taux,tauy,q,qr,qb,tpf(1,3)) elseif (initq .eq. 5) then c.....*new* Richard-Benno formulation using Solar Radiation & Clouds: call it_catch (nsst, tsst, nstep, it1, it2, tscl) if (it2 .ne. isst) then isst = it2 do i = 1, npt sstf(i,1) = sstf(i,2) cldf(i,1) = cldf(i,2) slrf(i,1) = slrf(i,2) enddo call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', tp, sstf(1,2)) call read_zt (idf_cld, lsst, npt, 1, it2, 'cld', tp, cldf(1,2)) call read_zt (idf_slr, lsst, npt, 1, it2, 'solr',tp, slrf(1,2)) endif do i = 1, npt sstf(i,3) = sstf(i,1) + tscl * (sstf(i,2) - sstf(i,1)) cldf(i,3) = cldf(i,1) + tscl * (cldf(i,2) - cldf(i,1)) slrf(i,3) = slrf(i,1) + tscl * (slrf(i,2) - slrf(i,1)) enddo call hflx_s94(npt,t,taux,tauy,sstf(1,3),cldf(1,3),slrf(1,3),q,qr,qb) elseif (initq .eq. 7) then tenso = enso_start + enso_scale * nstep call hflx_s89(tenso,npt,iox,t,sstf,cldf,ym,taux,tauy,q,qr,qb,tpf) elseif (initq .eq. 8) then c.....PBL model call it_catch (nsst, tsst, nstep, it1, it2, tscl) NXY = nx*ny newread = 0 if (it2 .ne. isst) then newread = 1 isst = it2 do i = 1, npt sstf(i,1) = sstf(i,2) cldf(i,1) = cldf(i,2) slrf(i,1) = slrf(i,2) wnsp(i,1) = wnsp(i,2) uwnd(i,1) = uwnd(i,2) vwnd(i,1) = vwnd(i,2) enddo do i = 1, NXY ahum(i,1) = ahum(i,2) atem(i,1) = atem(i,2) enddo call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', tp, sstf(1,2)) call read_zt (idf_cld, lsst, npt, 1, it2, 'cld', tp, cldf(1,2)) call read_zt (idf_slr, lsst, npt, 1, it2, 'solr',tp, slrf(1,2)) call read_zt (idf_wsp, lsst, npt, 1, it2, 'wndspd', tp, wnsp(1,2)) call read_zt (idf_uwd, lsst, npt, 1, it2, 'uwnd', tp, uwnd(1,2)) call read_zt (idf_vwd, lsst, npt, 1, it2, 'vwnd', tp, vwnd(1,2)) #ifdef gidon call read_zt (idf_ah, 0, NXY, 1, it2, 'spechum', tp, ahum(1,2)) call read_zt (idf_at, 0, NXY, 1, it2, 'Tair', tp, atem(1,2)) #else call read_zt (idf_ah, 0, NXY, 1, it2, 'airhum', tp, ahum(1,2)) call read_zt (idf_at, 0, NXY, 1, it2, 'airtem', tp, atem(1,2)) #endif endif if (newread.eq.1 .or.FIRST_STEP .or.mod(nstep, nstep_pbl).eq.0) then do i = 1, npt sstf(i,3) = sstf(i,1) + tscl * (sstf(i,2) - sstf(i,1)) cldf(i,3) = cldf(i,1) + tscl * (cldf(i,2) - cldf(i,1)) tpf(i,2) = uwnd(i,1) + tscl * (uwnd(i,2) - uwnd(i,1)) tpf(i,3) = vwnd(i,1) + tscl * (vwnd(i,2) - vwnd(i,1)) enddo do i = 1, npt wnd_speed = wnsp(i,1) + tscl * (wnsp(i,2) - wnsp(i,1)) if (wnd_speed .lt. pbl_wmin) then tpf(i,1) = pbl_wmin else tpf(i,1) = wnd_speed endif enddo do i = 1, NXY ahum(i,3) = ahum(i,1) + tscl * (ahum(i,2) - ahum(i,1)) atem(i,3) = atem(i,1) + tscl * (atem(i,2) - atem(i,1)) enddo call htflux_pbl (npt, nx, ny, iox, xm, ym, * t,cldf(1,3), tpf(1,1),tpf(1,2),tpf(1,3), ahum(1,3),atem(1,3), * qbf(1,2), qbf(1,3), qbf(1,4), amhum, amth) endif qcon_inv = 1./QCON qcon_gam = solr_gamma/QCON do i = 1, npt qsolr = slrf(i,1) + tscl * (slrf(i,2) - slrf(i,1)) qbf(i,1) = qsolr qbf(i,5) = 30.*(sstf(i,3) - t(i)) c...........Total Heat Flux at the surface: c...........Q = (1-gamma)*Q_sol - Q_lh - Q_sh - Q_lw qr(i) = qcon_gam * qsolr qtot = qsolr - qbf(i,2) - qbf(i,3) - qbf(i,4) q(i) = qcon_inv * qtot - qr(i) enddo endif return end c--------------------------------------------------------- subroutine epforc(nstep, npt, salt, sssf, evpf, qbf) c--------------------------------------------------------- c update EP using current sal(1), and SSS(i) implicit real(a-h,o-z),integer(i-n) include 'comm_data.h' include 'comm_new.h' dimension salt(1),sssf(npt,1),evpf(npt,1), qbf(npt,1) parameter (R_MMDAY2MSEC = 1./(24. * 3600. * 1000.)) parameter (CLATHT2EVAP = 1./(2.5e6*1028.)) if (initep .eq. 0) then do i = 1, npt evpf(i,1) = trans_coef * (SATM - salt(i)) enddo elseif (initep .eq. 1 .or. initep .eq. 2) then do i = 1, npt evpf(i,1) = trans_coef * (sssf(i,1) - salt(i)) enddo elseif (initep .eq. 3) then if (isss .ne. isst) then isss = isst do i = 1, npt sssf(i,1) = sssf(i,2) enddo call read_zt (idf_sss, lsst, npt, 1, isss, 'sss', tp, sssf(1,2)) endif do i = 1, npt sss_d = sssf(i,1) + sst_tscl * (sssf(i,2) - sssf(i,1)) evpf(i,1) = trans_coef * (sss_d - salt(i)) enddo elseif (initep .eq. 4) then call it_catch (nevp, tevp, nstep, it1, it2, tscl) if (it2 .ne. ievp) then ievp = it2 do i = 1, npt evpf(i,2) = evpf(i,3) enddo call read_zt (idf_evp, levp, npt, 1, it2, 'evp', tp, evpf(1,3)) endif do i = 1, npt c...........assumes that evp data in [mm/day]: E_P = evpf(i,2) + tscl * (evpf(i,3) - evpf(i,2)) evpf(i,1) = R_MMDAY2MSEC * salt(i) * E_P enddo elseif (initep.eq.8 .or. initep.eq.9) then call it_catch (nsst, tsst, nstep, it1, it2, tscl) if (it2 .ne. ievp) then ievp = it2 do i = 1, npt evpf(i,2) = evpf(i,3) enddo call read_zt (idf_evp, lsst, npt, 1, it2, 'precip', tp, evpf(1,3)) endif if (initep.eq.8) then do i = 1, npt precip = R_MMDAY2MSEC *(evpf(i,2) + tscl*(evpf(i,3) - evpf(i,2))) #ifdef gidon evapor = CLATHT2EVAP * qbf(i,2) * 0.7 #else evapor = CLATHT2EVAP * qbf(i,2) #endif evpf(i,1) = (evapor - precip) * salt(i) enddo else !...........see Gilles for the formula: do i = 1, npt precip = R_MMDAY2MSEC *(evpf(i,2) + tscl*(evpf(i,3) - evpf(i,2))) rlhvap = 2500800. - 2300. * t(i) evapor = qbf(i,2) / (rlhvap*(dens(i)+1000.)) evpf(i,1) = (evapor - precip) * salt(i) enddo endif endif return end c--------------------------------------------- subroutine hbcset(npt, nzp, nsig, lok, hmf, hclf) c--------------------------------------------- c apply B.C. to the H field if relaxing to climatology. implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' include 'comm_para.h' common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) dimension hmf(npt,1),hclf(npt,nzp,1) dimension lok(4*MAXSID,nz) nlo = nlok(1) if (icl_h .eq. 1) then do k = 1, nsig do n = 1, nlo i = lok(n,1) hmf(i,k) = hclf(i,k,1) enddo enddo elseif (icl_h .eq. 2) then do k = 1, nsig do n = 1, nlo i = lok(n,1) hmf(i,k) = hclf(i,k,1) + clm_tscl*(hclf(i,k,2) - hclf(i,k,1)) enddo enddo endif return end c----------------------------------------------------- subroutine tbcset(npt, nzp, lok, t0, hmf, tmf, tclf) c----------------------------------------------------- c apply B.C. to the temperature field. The land boundaries have zero c heat flux, the open ocean boundaries have a specified temperature. c c lok = (input) regular or compressed x-sort indices of the "open c ocean" boundary points at which t is constant. c nlok = (common) number of open ocean grid points. nlo .gt. 0 is c equivalent to having mtc=1 in previous versions c of the model. implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' include 'comm_para.h' common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) dimension t0(1),tmf(npt,1),hmf(npt,1),tclf(npt,nzp,1) dimension lok(4*MAXSID,nz) c if (icl_ts .eq. 0) then do k = 1, nz nlk = nlok(k) if (nlk.eq.0) return tik = t0(k) do n = 1, nlk i = lok(n,k) tmf(i,k) = tik * hmf(i,k) enddo enddo elseif (icl_ts .eq. 1) then do k = 1, nz nlk = nlok(k) if (nlk.eq.0) return do n = 1, nlk i = lok(n,k) tmf(i,k) = tclf(i,k,1)*hmf(i,k) enddo enddo elseif (icl_ts .eq. 2) then do k = 1, nz nlk = nlok(k) if (nlk.eq.0) return do n = 1, nlk i = lok(n,k) tmf(i,k) = hmf(i,k) * * (tclf(i,k,1) + clm_tscl*(tclf(i,k,2) - tclf(i,k,1))) enddo enddo endif return end c---------------------------------------------------------------------- subroutine it_catch (NN, tt, nstep, it1, it2, tscl) c---------------------------------------------------------------------- c.....Returns relative shift & indexes which are bracket nstep. implicit real(a-h,o-z),integer(i-n) dimension tt(1) include 'comm_new.h' denso = enso_start + enso_scale * nstep tstep = tt(2) - tt(1) if (tt(nn)-tt(1)+tstep .eq. 12.) then c.....Periodic Climatology Data denso = denso - 12.*(int(denso)/12) do it2 = 1, NN if (denso .lt. tt(it2)) goto 100 enddo denso = denso - 12. it2 = 1 100 if (it2 .eq. 1) then it1 = NN tscl = (12. - tt(NN) + denso)/(12. - tt(NN) + tt(1)) else it1 = it2 - 1 tscl = (denso - tt(it1))/(tt(it2) - tt(it1)) endif else c.....Non-periodic Data do it2 = 1, NN if (denso .lt. tt(it2)) goto 200 enddo it2 = NN+1 tscl = 0. 200 if (it2 .eq. 1) then it1 = 1 tscl = 0. elseif(it2 .eq. NN+1) then it1 = NN it2 = NN tscl = 0. else it1 = it2 - 1 tscl = (denso - tt(it1))/(tt(it2) - tt(it1)) endif endif return end c-------------------------------------------------------- subroutine data_on_model_grid (idf, lret, tag) c-------------------------------------------------------- character*(*) tag 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 include 'comm_new.h' include 'comm_data.h' logical grids_equiv lret = 1 if (idatgr .eq. 0) then c........check if data on the same grid as the model: if ( grids_equiv(idf, nxp,nyp,nxyc, nsx,nsy, tp)) lret = 0 c........check if data on the same grid as previous data: elseif ( grids_equiv(idf, mxp,myp,mpack,msx,msy, tp)) then lret = 1 c........check if data on the same grid as the model: elseif ( grids_equiv(idf, nxp,nyp,nxyc, nsx,nsy, tp)) then lret = 0 else write(6, *) tag, 'Only one data GRID allowed! Stop.' stop endif return end c------------------------------------------------------------------- logical function grids_equiv (idf, kxp,kyp,kpack,ksx,ksy, tpp) c------------------------------------------------------------------- dimension tpp(1), kask(1000) logical odb_ifatt common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc include 'comm_new.h' include 'comm_data.h' m_p = 0 call odb_rddm(idf, 'NPACK', m_p) m_x = 0 call odb_rdgr(idf, 'X', m_x, tpp) xer = 1.e-6*(xm(nxp) - xm(1)) if (xm(1) .lt. tpp(1)-xer .or. xm(nxp) .gt. tpp(m_x)+xer) then write (6, *) '!!! X grid of DATA must cover the model region' stop endif m_sx = 0 if ( odb_ifatt(idf, 'X', 'stretched') ) * call odb_getiattr(idf, 'X', 'stretched', m_sx) m_y = 0 yer = 1.e-6*(ym(nyp) - ym(1)) call odb_rdgr(idf, 'Y', m_y, tpp(m_x+1)) if (ym(1).lt. tpp(m_x+1)-yer .or. ym(nyp) .gt. tp(m_x+m_y)+yer) then write (6, *) '!!! Y grid of DATA must cover the model region' stop endif m_sy = 0 if ( odb_ifatt(idf, 'Y', 'stretched') ) * call odb_getiattr(idf, 'Y', 'stretched', m_sy) m_seg = 0 call odb_rddm (idf, 'NMASK', m_seg) call odb_rdvar(idf, 'MASK', tpp(m_x+m_y+1)) if (m_p.ne.kpack .or. * m_x.ne.kxp .or. m_y.ne.kyp .or. * m_sx.ne.ksx .or. m_sy.ne.ksy ) then grids_equiv = .FALSE. if (idatgr .eq. 0) then mpack = m_p mxp = m_x msx = m_sx myp = m_y msy = m_sy mseg = m_seg call datagrid_memory(tp) call blin_indx(tp) call blin_coef(tp(mseg+1)) endif else grids_equiv = .TRUE. endif return end c----------------------------------------------------------------- subroutine read_zt (idf, key, npt, iz, it, tag, ftmp, fdata) c----------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension ftmp(1), fdata(1) character*(*) tag include 'comm_data.h' if (key .eq. 0) then !! data on MODEL grid call odb_rd1v3 (idf, iz, it, tag, fdata) else !! data on a different grid call odb_rd1v3 (idf, iz, it, tag, ftmp) call blin_intr(npt, ixd, im2d, blcf, ftmp, fdata) endif return end c------------------------------------------------------- integer function nearest(ixy, xd, yd, iseg) c------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension xd(1), yd(1), iseg(1) include 'comm_new.h' iy0 = 1 + (ixy-1)/mxp ix0 = ixy - mxp*(iy0-1) x0 = xd(ix0) y0 = yd(iy0) nearest = iseg(1) iy = 1 + (nearest-1)/mxp ix = nearest - (iy-1)*mxp dmin = (xd(ix)-x0)**2 + (yd(iy)-y0)**2 do j = 2, mseg i = iseg(j) iy = 1 + (i-1)/mxp ix = i - mxp*(iy-1) d = (xd(ix)-x0)**2 + (yd(iy)-y0)**2 if (d .lt. dmin) then dmin = d nearest = i endif enddo return end c-------------------------------------- subroutine blin_coef(iw) c-------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension xd(1), yd(1), xm(1), ym(1), iw(1) common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc include 'comm_data.h' include 'comm_new.h' call bracket(mxp, xd, nxp, xm, iw) call bracket(myp, yd, nyp, ym, iw(nxp+1)) do k = 1, nxyc j = 1 + (iox(k)-1)/nxp i = iox(k) - (j-1)*nxp c.....find the i,j location for the four surrounding DATA grid points. i1 = iw(i) j1 = iw(nxp+j) c.....x-sort index of four DATA points surrounding MODEL point (i,j). im2d(k) = i1 + (j1-1)*mxp c.....find the interpolation ratios. fx = (xm(i)-xd(i1))/(xd(i1+1)-xd(i1)) fy = (ym(j)-yd(j1))/(yd(j1+1)-yd(j1)) blcf(k) = (1.-fx)*(1.-fy) blcf(k+nxyc) = fx*(1.-fy) blcf(k+2*nxyc) = (1.-fx)*fy blcf(k+3*nxyc) = fx*fy enddo return end c ----------------------------------------------------------------- subroutine bracket(nx1,x1,nx2,x2,it) c ----------------------------------------------------------------- c find the elements of x1 which bracket each element of x2. c returns it(i), for i=1,nx2 such that: c x1(it(i)) .le. x2(i) .and. x2(i) .le. x1(it(i)+1) c nx1 = (input) length of x1. c x1 = (input) must have x1(i+1) .gt. x1(i), i=1,nx1-1. c nx2 = (input) length of x2. c x2 = (input) must have x2(i+1) .gt. x2(i), i=1,nx2-1. c it = (output) nx2 indices of the lower side of the pair of c consecutive elements of x1 which bracket x2(i). c c must input x1(1).le.x2(1) .and. x1(nx1).ge.x2(nx2). c dimension x1(1),x2(1),it(1) c i1 = 1 do 20 i2=1,nx2 10 if(x2(i2).ge.x1(i1) .and. x2(i2).le.x1(i1+1)) goto 20 i1 = i1 + 1 if(i1.lt.nx1-1) goto 10 20 it(i2) = i1 return c end of bracket. end c---------------------------------------------------------------- subroutine clim_updt(npt,nz,nstep, h0,sigma,dzin,hclf,tclf,sclf,dclf) c---------------------------------------------------------------- include 'comm_new.h' include 'comm_data.h' dimension h0(1),hclf(npt,nz,1),tclf(npt,nz,1),sclf(npt,nz,1) dimension dclf(npt,nz),sigma(nz),dzin(nz+1) common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch if (icl_ts .ne. 2) return call it_catch (ntclm, tclm, nstep, it1, it2, clm_tscl) c.....if we still withis bracketing months: if (it2 .eq. iclm) return iclm = it2 if (icl_h .eq. 2 ) then !!..H-clim is time dependent do k = 1, nz !!...save old time do i = 1, npt hclf(i,k,1) = hclf(i,k,2) enddo enddo if (icl_htop .eq. 1) * call read_zt (idf_hcl, lclm,npt, 1, it2, 'mltc', tp, hclf(1,1,2)) sigk = sigma(3) do i = 1, npt hclf(i,2,2) = 0.5*hclf(i,1,2) + sigk*(z_begin - 1.5*hclf(i,1,2)) enddo do k = 3, nsig - 1 sigkp = sigma(k+1) do i = 1, npt hclf(i,k,2) = (sigk+sigkp) * (z_begin - 1.5*hclf(i,1,2)) enddo sigk = sigkp enddo k = nsig do i = 1, npt hclf(i,k,2) = dzin(k+1) + sigma(k)*(z_begin - 1.5*hclf(i,1,2)) enddo endif do k = 1, nz !!...save old time do i = 1, npt tclf(i,k,1) = tclf(i,k,2) enddo enddo if (use_salt) then do k = 1, nz !!...save old time do i = 1, npt sclf(i,k,1) = sclf(i,k,2) enddo enddo endif k2_h = 1 if (icl_h .eq. 2) k2_h = 2 call read_linz(idf_t, lclm, npt,mpack,nz,nzclm, it2, * hclf(1,1,k2_h),tclf(1,1,2),zclm,tp, 'temp') if (use_salt) then call read_linz(idf_s, lclm, npt,mpack,nz,nzclm, it2, * hclf(1,1,k2_h),sclf(1,1,2),zclm,tp, 'salt') endif if (ipre.eq.1) then if (use_salt) call potn_dens (npt,nzi,tclf(1,1,2),sclf(1,1,2),dclf) call dconv_cl (npt,nz,nzi,hclf,tclf(1,1,2),sclf(1,1,2),dclf) endif return end c----------------------------------------------------------------------- subroutine clim_relax(npt,nz,hmf,tmf,smf,hclf,tclf,sclf) c----------------------------------------------------------------------- include 'comm_new.h' include 'comm_data.h' dimension hmf(npt,1),tmf(npt,1),smf(npt,1), * hclf(npt,nz,1),tclf(npt,nz,1),sclf(npt,nz,1) if (icl_rlx .eq. 0) return !! NO RELAXATION - EXIT if (icl_h .eq. 1) then !! H-clim time independent do i = 1, npt coef = relax(i) do k = 1, nzi(i) hmf(i,k) = hmf(i,k) - coef*(hmf(i,k)-hclf(i,k,1)) enddo enddo elseif (icl_h .eq. 2) then !! h_clim time varying do i = 1, npt coef = relax(i) if (coef .ne. 0.) then do k = 1, nzi(i) hcl = hclf(i,k,1) + clm_tscl*(hclf(i,k,2) - hclf(i,k,1)) hmf(i,k) = hmf(i,k) - coef*(hmf(i,k) - hcl) enddo endif enddo endif if (icl_ts .eq. 1) then if (use_salt) then do i = 1, npt coef = relax(i) if (coef .ne. 0.) then do k = 1, nzi(i) tmf(i,k) = tmf(i,k) - coef*(tmf(i,k)-tclf(i,k,1)) smf(i,k) = smf(i,k) - coef*(smf(i,k)-sclf(i,k,1)) enddo endif enddo else do i = 1, npt coef = relax(i) if (coef .ne. 0.) then do k = 1, nzi(i) tmf(i,k) = tmf(i,k) - coef*(tmf(i,k)-tclf(i,k,1)) enddo endif enddo endif elseif (icl_ts .eq. 2) then !!...vary hmix/hthermo if (use_salt) then do i = 1, npt coef = relax(i) if (coef .ne. 0.) then do k = 1, nzi(i) tcl = tclf(i,k,1) + clm_tscl*(tclf(i,k,2) - tclf(i,k,1)) scl = sclf(i,k,1) + clm_tscl*(sclf(i,k,2) - sclf(i,k,1)) tmf(i,k) = tmf(i,k) - coef*(tmf(i,k) - tcl) smf(i,k) = smf(i,k) - coef*(smf(i,k) - scl) enddo endif enddo else do i = 1, npt coef = relax(i) if (coef .ne. 0.) then do k = 1, nzi(i) tcl = tclf(i,k,1) + clm_tscl*(tclf(i,k,2) - tclf(i,k,1)) tmf(i,k) = tmf(i,k) - coef*(tmf(i,k) - tcl) enddo endif enddo endif endif return end c--------------------------------- subroutine afill(n, a, v) c--------------------------------- dimension a(1) do i = 1, n a(i) = v enddo return end c----------------------------------------------------------------------------- subroutine read_linz(idf,key,NPT,MPT,NZ,MZ,it, hdat,fdat,zvert,fvert,tag) c----------------------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension hdat(npt,1),fdat(npt,1), zvert(1),fvert(1) character*(*) tag include 'comm_new.h' include 'comm_data.h' dimension aa(npt,1), bb(mpt,1) pointer (p_aa, aa), (p_bb, bb) if (key .eq. 0) then call mem_alloc(p_aa, MZ*npt, 2, 'AA space in read_linz') do k = 1, MZ call odb_rd1v3(idf, k, it, tag, aa(1,k)) enddo do i = 1, npt do k = 1, mz fvert(k) = aa(i,k) enddo call zlin_intrp (i, npt,nz,mz, hdat,fdat,zvert,fvert) enddo call mem_free(p_aa, MZ*npt, 2) else call mem_alloc(p_bb, MZ*mpt, 2, 'BB space in read_linz') do k = 1, MZ call odb_rd1v3(idf, k, it, tag, bb(1,k)) enddo call zlin_blin(NPT,MPT,NZ,MZ,ixd,im2d,blcf,bb,hdat,fdat,zvert,fvert) call mem_free(p_bb, MZ*mpt, 2) endif return end c----------------------------------------------------------------------- subroutine zlin_blin(NPT,MPT,NZ,MZ,ixd,im2d,blcf,aa,h,f,zval,fval) c----------------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' dimension ixd(1),im2d(1),blcf(npt,1) dimension aa(mpt,1), h(npt,1),f(npt,1), zval(1),fval(1) do i = 1, npt i1 = im2d(i) n1 = ixd(i1) n2 = ixd(i1+1) n3 = ixd(i1+mxp) n4 = ixd(i1+mxp+1) b1 = blcf(i,1) b2 = blcf(i,2) b3 = blcf(i,3) b4 = blcf(i,4) do k = 1, mz fval(k) = b1*aa(n1,k) + b2*aa(n2,k) + b3*aa(n3,k) + b4*aa(n4,k) enddo call zlin_intrp (i, npt,nz,mz, h,f,zval,fval) enddo return end c------------------------------------------------------------------------ subroutine blin_intr(npt, ixd, im2d, blcf, fd, f) c------------------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) dimension ixd(1), im2d(1), blcf(npt,1), f(1), fd(1) include 'comm_new.h' do i = 1, npt i1 = im2d(i) f(i) = blcf(i,1)*fd(ixd(i1)) + blcf(i,2)*fd(ixd(i1+1)) * + blcf(i,3)*fd(ixd(i1+mxp)) + blcf(i,4)*fd(ixd(i1+mxp+1)) enddo return end c------------------------------------------------------- subroutine blin_indx (iseg) c------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension iseg(2,1) include 'comm_data.h' include 'comm_new.h' c.....fill in ixd() according with the data compression iseg(2,mseg/2): k = 0 do i = 1, mseg/2 do j = iseg(1,i), iseg(2,i) k = k + 1 ixd(j) = k enddo enddo c.....continue data to all points do k = 1, mxp*myp if (ixd(k) .eq. 0) ixd(k) = ixd(nearest(k, xd, yd, iseg)) enddo return end c------------------------------------------------------------- subroutine zlin_intrp(i, npt, nz, mz, h, f, zval, fval) c------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension h(npt,1), f(npt,1), zval(1), fval(1) zbot = zval(mz) fbot = fval(mz) mnext = 2 scale = (fval(mnext)-fval(mnext-1))/(zval(mnext)-zval(mnext-1)) shift = fval(mnext-1) - scale*zval(mnext-1) dlay = 0.5*h(i,1) zlay = dlay do k = 1, nz if (zlay .gt. zval(mnext)) then if (zlay .gt. zbot) then shift = fbot scale = 0. else do m = mnext+1, mz if (zlay .le. zval(m)) then scale = (fval(m) - fval(m-1))/(zval(m) - zval(m-1)) shift = fval(m-1) - scale*zval(m-1) mnext = m goto 100 endif enddo endif endif 100 f(i,k) = shift + scale * zlay dlay = h(i,k) - dlay zlay = zlay + 2.*dlay enddo return end