subroutine setophys c c Sets up ocean circulation model. This is a combination of old c setogrd.F and initoce.F, without unused options and without c reading the restart files, and finally without mapping, which is c done elsewhere. c c#include "gloparam.F.h" #include "oparams.F.h" #include "ogrid.F.h" #include "ostate.F.h" #include "mapping.h" #include "radpth.h" real Wk1(len) c c==> Grid data file info real xse( IM+1,JM+1) ! Long of se corner of grid box (degrees) real yse( IM+1,JM+1) ! Lat of se corner of grid box (degrees) c real chi( IM,JM ) ! angle between coord and zonal dir (Radians) real Um ( IM,JM ) ! U metric term real Vm ( IM,JM ) ! V metric term real zx0 ! Location of center of restart grid box (1,1) real zy0 ! Location of center of restart grid box (1,1) real zdx ! Separation of restart grid boxes real zdy ! Separation of restart grid boxes real xc(IM,JM) real yc(IM,JM) real dxg( IM, JM+1) real dyg( IM+1,JM ) real dxgu( IM, JM) real dygu( IM,JM ) real cos_chig( IM, JM) real sin_chig( IM,JM ) real ag( IM, JM) real agu( IM, JM) real dx1,dx2,dx3 real tag c==> Functions real SSUM real G_CIRCLE_DIST integer ISMIN, ISMAX c c--------------------------------------------------------------------- print *,'#################################################' print *,'# Start SETOPHYS: ' c zx0 = along0 + 0.5*dlong zy0 = alat0 + 0.5*dlat zdx = dlong zdy = dlat write(6,*)'zx0,zy0,zdx,zdy = ',zx0,zy0,zdx,zdy do j=1,JM+1 do i=1,IM+1 xse(i,j) = zx0 + ( i - 1.5 ) * zdx yse(i,j) = zy0 + ( j - 1.5 ) * zdy enddo enddo write(6,*)'xse1,xse2 = ',xse(1,1),xse(2,1) write(6,*)'yse1,yse2 = ',yse(1,1),yse(1,2) do j=1,JM do i=1,IM c chi(i,j) = 0. Um(i,j) = TAN( yse(i+1,j+1) * DTOR ) / RADIUS Vm(i,j) = 0. enddo enddo print *,'# setting up ocean grid' ! The input grid has lat and long of corners defined ! This is IM+1 by JM+1 ! Get the lengths of each side ! Get the center long and lats ! Get the areas of each box ! Get the distances between centers ! Find grid box centers from corners ! Watch out for 360 degree longitude shifts do j=1,JM do i=1,IM dx1 = xse(i+1,j ) - xse(i,j) dx3 = xse(i ,j+1) - xse(i,j) dx2 = xse(i+1,j+1) - xse(i,j) dx1 = AMOD( dx1 + 540. , 360. ) - 180. dx2 = AMOD( dx2 + 540. , 360. ) - 180. dx2 = AMOD( dx3 + 540. , 360. ) - 180. xc(i,j) = xse(i,j) + ( dx1 + dx2 + dx3 )*.25 yc(i,j) = (yse(i,j) + yse(i+1,j) + yse(i,j+1) + yse(i+1,j+1))*.25 enddo enddo do j=1,JM do i=1,IM c cos_chig(i,j) = COS( chi(i,j) ) c sin_chig(i,j) = SIN( chi(i,j) ) cos_chig(i,j) = 1.0 sin_chig(i,j) = 0.0 enddo enddo ! Length of South box edge do j=1,JM+1 do i=1,IM dxg(i,j) = G_CIRCLE_DIST( xse(i ,j)*DTOR, yse(i ,j)*DTOR, & xse(i+1,j)*DTOR, yse(i+1,j)*DTOR, & RADIUS ) enddo enddo i = ISMIN( IM*(JM+1), dxg, 1 ) j = ISMAX(IM*(JM+1), dxg, 1 ) print *,'# Min, max of dxg = ',dxg(i,1),dxg(j,1) ! Length of west box edge do j=1,JM do i=1,IM+1 dyg(i,j) = G_CIRCLE_DIST( xse(i,j )*DTOR, yse(i,j )*DTOR, & xse(i,j+1)*DTOR, yse(i,j+1)*DTOR, & RADIUS ) enddo enddo i = ISMIN( (IM+1)*JM, dyg, 1 ) j = ISMAX( (IM+1)*JM, dyg, 1 ) print *,'# Min, max of dyg = ',dyg(i,1),dyg(j,1) ! EW dist between box centers do j=1,JM do i=1,IM-1 dxgu(i,j) = G_CIRCLE_DIST( xc(i ,j)*DTOR, yc(i ,j)*DTOR, & xc(i+1,j)*DTOR, yc(i+1,j)*DTOR, & RADIUS ) enddo enddo do j=1,JM dxgu(IM,j) = G_CIRCLE_DIST( xc(IM ,j)*DTOR, yc(IM ,j)*DTOR, & xc(1 ,j)*DTOR, yc(1 ,j)*DTOR, & RADIUS ) enddo ! NS dist between box centers do j=1,JM-1 do i=1,IM dygu(i,j) = G_CIRCLE_DIST( xc(i,j )*DTOR, yc(i,j )*DTOR, & xc(i,j+1)*DTOR, yc(i,j+1)*DTOR, & RADIUS ) enddo enddo ! Area of grid boxes do j=1,JM do i=1,IM ag(i,j) = 0.25*(dxg(i,j) + dxg(i,j+1)) * (dyg(i,j)+dyg(i+1,j)) enddo enddo i = ISMIN( IM*JM, ag, 1 ) j = ISMAX( IM*JM, ag, 1 ) print *,'# Min, max of ag = ',ag(i,1),ag(j,1) ! Area of u grid boxes do j=1,JM-1 do i=1,IM ip = MOD( i, IM ) + 1 agu(i,j) = 0.25*(ag(i,j) + ag(i,j+1) + ag(ip,j) + ag(ip,j+1)) enddo enddo i = ISMIN( IM*(JM-1), agu, 1 ) j = ISMAX( IM*(JM-1), agu, 1 ) print *,'# Min, max of agu = ',agu(i,1),agu(j,1) l = nwater+nland do i=1,l if (ieast(i).lt.1 .or. ieast(i).gt.l) then print *,'# Error in output from MAPPING, ieast failure:' print *,'# ',i,ieast(i),' E' stop endif if (iwest(i).lt.1 .or. iwest(i).gt.l) then print *,'# Error in output from MAPPING, iwest failure:' print *,'# ',i,iwest(i),' W' stop endif if (inorth(i).lt.1 .or. inorth(i).gt.l) then print *,'# Error in output from MAPPING, inorth failure:' print *,'# ',i,inorth(i),' N' stop endif if (isouth(i).lt.1 .or. isouth(i).gt.l) then print *,'# Error in output from MAPPING, isouth failure:' print *,'# ',i,isouth(i),' S' stop endif enddo do i=1,LEN Bh(i) = 0 Bu(i) = 0 Area_h(i) = 0. Area_h_r(i) = 0.0 Area_u(i) = 0. Area_u_r(i) = 0.0 oc_lat(i) = -999. oc_lon(i) = -999. oc_lat_ne(i) = -999. oc_lon_ne(i) = -999. Dx_h(i) = 0. Dy_h(i) = 0. U_metric(i) = 0. V_metric(i) = 0. Coriolis(i) = 0. Zbottom(i) = 0.0 cos_chi(i) = 0. sin_chi(i) = 0. enddo do i=1,nwater Bh(i) = 1 enddo do i=1,nwater Bu(i) = Bh(i)*Bh(ieast(i))*Bh(inorth(i))*Bh(ieast(inorth(i))) oc_lon(i) = xc ( gridi(i)+1, gridj(i)+1 ) oc_lat(i) = yc ( gridi(i)+1, gridj(i)+1 ) oc_lon_ne(i) = xse( gridi(i)+2, gridj(i)+2 ) oc_lat_ne(i) = yse( gridi(i)+2, gridj(i)+2 ) Area_h(i) = ag ( gridi(i)+1, gridj(i)+1 ) * Bh(i) Dx_h(i) = dxg( gridi(i)+1, gridj(i)+2 ) Dy_h(i) = dyg( gridi(i)+2, gridj(i)+1 ) Area_u(i) = agu( gridi(i)+1, gridj(i)+1 ) * Bu(i) U_metric(i) = um ( gridi(i)+1, gridj(i)+1 ) * Bu(i) V_metric(i) = vm ( gridi(i)+1, gridj(i)+1 ) * Bu(i) Coriolis(i) = 2.*OMEGA*SIN( oc_lat_ne(i)*DTOR ) * Bu(i) Zbottom(i) = Depth( gridi(i)+1, gridj(i)+1 ) * Bh(i) cos_chi(i) = cos_chig( gridi(i)+1, gridj(i)+1 ) * Bu(i) sin_chi(i) = sin_chig( gridi(i)+1, gridj(i)+1 ) * Bu(i) Area_h_r(i) = Bh(i) / ( Area_h(i) + 1.E-36) Area_u_r(i) = Bu(i) / ( Area_u(i) + 1.E-36) enddo totalah = SSUM(nwater,Area_h,1) tag = SSUM(IM*JM , Ag, 1) ! Printout grid statistics print *,'# POSEIDON GRID HAS BEEN SET UP FROM RESTART FILE' print *,'#' print 1, tag print 2, totalah print 3, 100.*totalah/tag print 4, 100.*totalah/( 4 * PI * RADIUS * RADIUS ) print *,'#' 1 format(' # The total area on the rectangular grid is ',1PE20.12, & ' sq. meters') 2 format(' # The total area in water is ',1PE20.12) 3 format(' # Or ',F8.3,' percent of the grid area') 4 format(' # Or ',F8.3,' percent of Earth surface') print *,'# done with ocean grid setup' print *,'#################################################' c Now adjust D_min and D_max to accomodate the bottom topography if ( EXTMODE ) then write (6,*) '# USING EXTERNAL MODE VERSION' do i=1,nwater D_max(i,km) = -Zbottom(i) ! Max depth (positive down) D_min(i,km) = -Zbottom(i) ! Z is positive UPwards enddo do k=km-1 ,1, -1 do i=1,nwater D_max(i,k) = AMIN1(D_max1d(k), D_max(i,k+1) - hmin_b(k)) D_min(i,k) = AMIN1(D_min1d(k), D_min(i,k+1) - hmin_b(k)) enddo enddo do k=1,km do i=nwater+1,LEN D_max(i,k) = 0. D_min(i,k) = 0. enddo enddo else ! EXTMODE do k=1,km do i=1,nwater D_min(i,k) = D_min1d(k) D_max(i,k) = D_max1d(k) enddo do i=nwater+1,len D_min(i,k) = 0. D_max(i,k) = 0. enddo enddo ! k endif ! EXTMODE c c do k=1,km c i=ISMIN( nwater, D_min(1,k), 1) c j=ISMAX( nwater, D_min(1,k), 1) c write (6,*) '# Min, Max of D_min(',k,') =',D_min(i,k),D_min(j,k) c i=ISMIN( nwater, D_max(1,k), 1) c j=ISMAX( nwater, D_max(1,k), 1) c write (6,*) '# Min, Max of D_max(',k,') =',D_max(i,k),D_max(j,k) c do i=1,nwater c Wk1(i) = D_max(i,k) + Zbottom(i) c enddo c i=ISMIN( nwater, Wk1, 1) c j=ISMAX( nwater, Wk1, 1) c write (6,*) '# Min, Max below D_max(',k,') =',Wk1(i),Wk1(j) c enddo c c Set up initial buoyancy do k=1,KM do i=1,nwater B(i,k,1) = BYNCY( T(i,k,1), S(i,k,1), 0.) B(i,k,2) = BYNCY( T(i,k,2), S(i,k,2), 0.) enddo do i=nwater+1,nwater+nland B(i,k,1) = 0. B(i,k,2) = 0. enddo enddo call PRNTAVGS( LEN, KM, nwater, extmode, & H, U, V, T, S, B, Area_h, Zbottom ) c c Initialize entrainment velocity do k = 1,km do i = 1,len W_ent(i,k) = 0.0 enddo enddo c print *,' done with ocean initialization' return end