C ********************************************************************** C * * C * SOFTWARE LICENSING * C * * C * This program is free software; you can redistribute * C * it and/or modify it under the terms of the GNU * C * General Public License as published by the Free * C * Software Foundation, either Version 2 of the * C * license, or (at your option) any later version. * C * * C * This program is distributed in the hope that it * C * will be useful, but without any warranty; without * C * even the implied warranty of merchantability or * C * fitness for a particular purpose. See the GNU * C * General Public License for more details. * C * * C * A copy of the GNU General Public License is * C * available at http://www.gnu.org/copyleft/gpl.html * C * or by writing to the Free Software Foundation, Inc.,* C * 59 Temple Place - Suite 330, Boston, MA 02111, USA. * C * * C ********************************************************************** SUBROUTINE SETDOM C C VERSION(06/27/91) C INCLUDE 'comdeck' C C FOR WET GRID INPUT FROM POM C CHARACTER*10 RESTAR COMMON /COAST1/ICNT,INDX(MAXWET),JNDX(MAXWET),RESTAR C SAVE C C----------------------------------------------------------------------- C INITIAL CONDITION PREPARATION PROGRAM FOR A GENERAL CIRCULATION C MODEL IN ORTHOGONAL CURVILINEAR COORDINATES C----------------------------------------------------------------------- C c DIMENSION ADVUA(IM,JM),ADVVA(IM,JM),ADVUU(IM,JM),ADVVV(IM,JM) c DIMENSION DRHOX(IM,JM,KB),DRHOY(IM,JM,KB),TRNU(IM,JM),TRNV(IM,JM) DIMENSION IVAR(IM,JM),PRT(IM,KB) EQUIVALENCE (IVAR,TPS),(PRT,A) DIMENSION COM(80) DIMENSION ISTDAM(NUMDAM),JSTDAM(NUMDAM),NTDAM(NUMDAM) CHARACTER DIRTDAM(NUMDAM)*4 C C-------- ALL UNITS IN M.K.S. SYSTEM ----------------------------------- ISTART=1 IEND=1 INTX=0 C C FOR EXTERNAL HYDRODYNAMICS RUN C C ASSUME POM OUTPUT FOR WET GRID C IF (HYDTYPE.EQ.'EXTERNAL') THEN IF(RESTAR.EQ.'HOT START '.AND.IWET.EQ.1) GOTO 153 cqa OPEN (IUTRN,FILE='gcm_geom',FORM='UNFORMATTED') READ (IUTRN)DZ,DZZ READ (IUTRN)H,H1,H2,TPS READ (IUTRN) ANG,NU C CONVERTED TO DEGREES FROM RADIANS BY QUAMRUL QA 9/15/98 DO J=1,JM DO I=1,IM ANG(I,J)=ANG(I,J)*360.0/(2.*3.141593) ENDDO ENDDO CLOSE(IUTRN) C OPEN (IUTRN,FORM='unformatted',FILE='wet_grid') READ (IUTRN)ICNT DO I=1,ICNT READ(IUTRN)INDX(I),JNDX(I) END DO CLOSE (IUTRN) C DO 151 K=1,KBM1 DZR(K)=1./DZ(K) 151 CONTINUE C C RESET LAND ELEMENTS TO ZERO DEPTH (SET TO 1.000 m IN POM) C DO 152 J=1,JM DO 152 I=1,IM IF (H(I,J).LT.1.002) H(I,J)=0.0 152 CONTINUE GOTO 999 C cqa cqa Reads wet_grid when "HOT START" and IWET=1 i.e. IWET on cqa 153 CONTINUE OPEN (IUTRN,FORM='unformatted',FILE='wet_grid') READ (IUTRN)ICNT DO I=1,ICNT READ(IUTRN)INDX(I),JNDX(I) END DO CLOSE (IUTRN) RETURN ENDIF C OPEN(IUGRD,FILE='model_grid') C WRITE(IUPRT,51) 51 FORMAT(/'...... MODEL STARTING UP FROM INITAL CONDITIONS .......') C C-------- ESTABLISH DEPTH ARRAY ---------------------------------------- READ(IUGRD,11) (COM(I),I=1,80) WRITE(IUPRT,3) (COM(I),I=1,80) 3 FORMAT(/1X,80A1/) 11 FORMAT(80A1) C C-------- INITIALIZE SIGMA LEVELS -------------------------------------- READ(IUGRD,11) (COM(I),I=1,80) WRITE(IUPRT,3) (COM(I),I=1,80) C READ(IUGRD,4) IKB 4 FORMAT(I5) WRITE(IUPRT,40) IKB 40 FORMAT(' KB = ',I5,/) IF(IKB.GT.101) THEN WRITE(6,41) KBM1 CALL SYSTEM ('rm gcm_temp*') STOP ENDIF 41 FORMAT(//' NUMBER OF MODEL LAYERS',I5,' (KBM1)'/ . ' GREATER THAN'/ . ' MAXIMUM NUMBER OF VERTICAL DISTRIBUTION LAYERS'/ . ' SPECIFIED IN run_data DISCHARGE INFORMATION'/ . ' PLEASE CORRECT THIS PROBLEM AND TRY AGAIN'//) C IF(IKB.NE.KB) THEN WRITE(6,42) IKB,KB CALL SYSTEM ('rm gcm_temp*') STOP ENDIF 42 FORMAT(//' NUMBER OF SIGMA LEVELS IN MODEL_GRID',I5,' (IKB)'/ . ' NOT EQUAL TO'/ . ' NUMBER OF SIGMA LEVELS IN COMDECK ',I5,' (KB)'/ . ' PLEASE CORRECT THIS PROBLEM AND TRY AGAIN'//) C IF(IKB.LT.3 .AND. TOR.NE.'BAROTROPIC') THEN WRITE(6,43) IKB,TOR CALL SYSTEM ('rm gcm_temp*') STOP ENDIF 43 FORMAT(//' IKB/KB = ',I5,' NOT ALLOWED IN A ',A10,' RUN'/ . ' NUMBER OF SIGMA LEVELS MUST BE GREATER THAN 2'/ . ' PLEASE CORRECT THIS PROBLEM AND TRY AGAIN'//) C DO 140 K=1,IKB READ(IUGRD,5) Z(K) 140 CONTINUE WRITE(IUPRT,5) (Z(K),K=1,IKB) 5 FORMAT(8F10.5) C DO 150 K=1,KBM1 DZ(K)=Z(K)-Z(K+1) DZR(K)=1./DZ(K) ZZ(K)=.5*(Z(K)+Z(K+1)) 150 CONTINUE DO 460 K=1,KBM2 DZZ(K)=ZZ(K)-ZZ(K+1) 460 CONTINUE DZZ(KBM1)=0.0 DZ(KB)=0.0 C C-------- DEFINE THE METRICS OF THE COORDINATE TRANSFORMATION ---------- READ(IUGRD,11) (COM(I),I=1,80) WRITE(IUPRT,3) (COM(I),I=1,80) C READ(IUGRD,7) IIX, IJY WRITE(IUPRT,71) IIX, IJY 7 FORMAT(2I5,6F10.2) 91 FORMAT(2I5,4F10.2,2F10.5,f5.1) c 91 FORMAT(2I4,8f9.0) 71 FORMAT(' IM = ',I5,/' JM = ',I5) IF(IIX.NE.IM) THEN WRITE(6,8) IIX,IM CALL SYSTEM ('rm gcm_temp*') STOP ENDIF 8 FORMAT (//' MODEL_GRID I-INDEX',I5,' (IIX)',/ . ' DOES NOT EQUAL'/ . ' COMDECK I-INDEX',I5,' (IM)'/ . ' PLEASE CORRECT THIS PROBLEM AND TRY AGAIN'//) IF(IJY.NE.JM) THEN WRITE(6,9) IJY,JM CALL SYSTEM ('rm gcm_temp*') STOP ENDIF 9 FORMAT (//' MODEL_GRID J-INDEX',I5,' (IJY)',/ . ' DOES NOT EQUAL'/ . ' COMDECK J-INDEX',I5,' (JM)'/ . ' PLEASE CORRECT THIS PROBLEM AND TRY AGAIN'//) C C..........SET PARAMETERS TO INITIAL VALUES.......... DO 50 J=1,JM DO 50 I=1,IM H1(I,J) = 0.0 H2(I,J) = 0.0 H(I,J) = 0.0 ANG(I,J) = 0.0 COR(I,J) = 0.0 DATUM(I,J) = 0.0 YGRID(I,J) = 0.0 50 CONTINUE C NUMB=IIX*IJY DO 100 N=1,NUMB cqa 100 READ(IUGRD,7,END=101) cqa . I,J,H1(I,J),H2(I,J),H(I,J),ANG(I,J),COR(I,J),DATUM(I,J) READ(IUGRD,91,END=102,ERR=101) ! READ(IUGRD,92,END=102,ERR=101) !RMarsooli_Jan2015 . I,J,H1(I,J),H2(I,J),H(I,J),ANG(I,J),COR(I,J),XGRID(I,J), . DATUM(I,J) YGRID(I,J)=COR(I,J) CNG IF(J.LE.8) H(I,J)=17. ! BIG HARDWIRE HERE H(I,J)=H(I,J)+WETMIN ! CNG09082011 MAKE GRID DEEPER BY WETMIN 100 continue ! 92 FORMAT(2I5,4F10.2,2F10.5,f5.3) !RMarsooli_Jan2015 101 CONTINUE BACKSPACE(IUGRD) READ(IUGRD,11,END=102) (COM(I),I=1,80) WRITE(IUPRT,3) (COM(I),I=1,80) READ(IUGRD,7)NUMTDAM IF(NUMTDAM.GT.NUMDAM)THEN WRITE(IUPRT,*)'NUMBER OF THIN DAM IS GREATER THAN ',NUMDAM WRITE(IUPRT,*)'PLEASE INCREASE THE ARRAY SIZE IN comdeck' STOP END IF READ(IUGRD,49) .(ISTDAM(N),JSTDAM(N),DIRTDAM(N),NTDAM(N),N=1,NUMTDAM) 49 FORMAT(2I5,1X,A4,I5) 102 CONTINUE CLOSE(IUGRD) call RDCORLOC !csvv080327 C C-------- DEFINE MASK FOR FREE SURFACE HEIGHT = FSM -------------------- C-------- DEFINE MASK FOR (U,V) VELOCITY = (DUM,DVM) ------------------- C-------- IF DEPTH <= 0.0, DEPTH AND MASK SET TO 0 --------------------- 999 DO 400 J=1,JM DO 400 I=1,IM FSM(I,J)=1.0 DUM(I,J)=1.0 DVM(I,J)=1.0 IF(H1(I,J).LE.0.0) THEN ! BASE IT ON H1, NOT H, W&D NG04212010 ! IF(H(I,J).LE.0.0) THEN ! BASE IT ON H1, NOT H, W&D NG04212010 CPL 8/16/99 H(I,J)=0.0 cpl H(I,J)=0.0+1.0E-3 H(I,J)=0.0+1.0E-10 FSM(I,J)=0.0 DUM(I,J)=0.0 DVM(I,J)=0.0 ! ELSEIF(H(I,J).EQ.-999.99) THEN ! FLAG FOR OPEN BC CORNERS ! H(I,J)=0.0+1.0E-10 ! FSM(I,J)=0.0 ! DUM(I,J)=0.0 ! DVM(I,J)=0.0 ELSEIF(H(I,J).EQ.0.0) THEN H(I,J)=0.0+1.0E-10 ! ELSEIF(H(I,J).LT.-20.0) THEN ! H(I,J)=0.0+1.0E-10-1.*float(I*J)/float(IM*JM) ENDIF 400 CONTINUE CPL 8/16/99 400 H(I,J)=H(I,J)+1.0E-03 C ! DO 450 I=1,IM ! DUM(I,3)=0.0 !!!!!!!!!!!!!!!!!!!!!!!!!!! BIG HARDWIRE HERE ! 450 CONTINUE DO 420 J=1,JMM1 DO 420 I=1,IM IF(FSM(I,J).EQ.0.0.AND.FSM(I,J+1).NE.0.0) DVM(I,J+1)=0.0 420 CONTINUE DO 440 J=1,JM DO 440 I=1,IMM1 IF(FSM(I,J).EQ.0.0.AND.FSM(I+1,J).NE.0.0) DUM(I+1,J)=0.0 440 CONTINUE C FOR THIN DAM DO 180 N=1,NUMTDAM IS=ISTDAM(N) JS=JSTDAM(N) IF(DIRTDAM(N).EQ.'IDIR') THEN DO 181 I=IS,IS+NTDAM(N)-1 DVM(I,JS)=0.0 181 CONTINUE ELSE ! JDIR DO 182 J=JS,JS+NTDAM(N)-1 DUM(IS,J)=0.0 182 CONTINUE END IF 180 CONTINUE C DO 160 J=1,JM DO 160 I=1,IM DHMWAD(I,J)=FSM(I,J) !RMarsooli_June2015 DUMWAD(I,J)=DUM(I,J) ! NG05192010 DVMWAD(I,J)=DVM(I,J) ! NG05192010 H1(I,J)=H1(I,J)+1.E-10 H2(I,J)=H2(I,J)+1.E-10 ANG(I,J)=ANG(I,J)*2.*3.141593/360. COR(I,J)=2.*7.292E-5*SIN(COR(I,J)*2.*3.141593/360.)*FSM(I,J) CNG COR(I,J)=0.0 ! HARDWIRE COSANG(I,J)=COS(ANG(I,J)) !RMarsooli SINANG(I,J)=SIN(ANG(I,J)) !RMarsooli DELX(I,J)=SQRT((H1(I,J)*COSANG(I,J))**2+(H2(I,J)*SINANG(I,J))**2) !RMarsooli DELY(I,J)=SQRT((H1(I,J)*SINANG(I,J))**2+(H2(I,J)*COSANG(I,J))**2) !RMarsooli crm write(130,'(5f15.7)') ANG(I,J),H1(I,J),DELX(I,J),H2(I,J),DELY(I,J) 160 CONTINUE C !FLOODIT=AMAX1(AMIN1(WETEPS,WETMIN/2.),1.E-3) ! W&D NG04212010 MAKE SURE D>0 CNG FLOODIT=AMAX1(WETMIN-1.E-5,1.E-3) ! W&D NG04212010 MAKE SURE D>0 FLOODIT=AMAX1(WETMIN-1.E-3,1.E-3) ! W&D NG04212010 MAKE SURE D>0 DO 260 J=1,JM DO 260 I=1,IM IF(FSM(I,J).GT.0.0) THEN ! W&D NG05182010 IF ORIGINALLY ABOVE GEOID, OR JUST BELOW, MAKE SURE D>0=FLOODIT D(I,J)=AMAX1(FLOODIT,H(I,J)) DT(I,J)=D(I,J) ELB(I,J)=D(I,J)-H(I,J) EL(I,J)=D(I,J)-H(I,J) ETB(I,J)=D(I,J)-H(I,J) ET(I,J)=D(I,J)-H(I,J) EGB(I,J)=D(I,J)-H(I,J) CNG ELF(I,J)=D(I,J)-H(I,J) CNG ETF(I,J)=D(I,J)-H(I,J) IF (ELB(I,J).GT.0.0) WRITE(IUPRT,170) I,J,D(I,J) ELSE ! NG, JUST FOR TIME STEP CALCULATION D(I,J)=H(I,J) DT(I,J)=D(I,J) ENDIF 260 CONTINUE 170 FORMAT(' IJ : ',2I5,' BELOW GEOID. INITIAL D SET TO ',F7.3) C DO 340 J=2,JMM1 DO 340 I=2,IMM1 ART(I,J)=H1(I,J)*H2(I,J) ARU(I,J)=.25E0*(H1(I,J)+H1(I-1,J))*(H2(I,J)+H2(I-1,J)) ARV(I,J)=.25E0*(H1(I,J)+H1(I,J-1))*(H2(I,J)+H2(I,J-1)) 340 CONTINUE DO 310 J=1,JM ARU(IM,J)=ARU(IMM1,J) 310 ARU(1,J)=ARU(2,J) DO 320 I=1,IM ARV(I,JM)=ARV(I,JMM1) 320 ARV(I,1)=ARV(I,2) C DO 410 J=1,JM DO 410 I=1,IM 410 TPS(I,J)=SQRT(H1(I,J)**2+H2(I,J)**2) C CALL MAXMIN(TPS,FSM,IM,JM,FMAX,FMIN) C IF(TOR.NE.'BAROTROPIC') THEN DO 350 K=1,KBM1 DO 350 J=1,JM DO 350 I=1,IM AAM(I,J,K)=HORCON*VARHF(I,J)*TPS(I,J)/FMIN*FSM(I,J) 350 CONTINUE DO 375 J=1,JM DO 375 I=1,IM 375 AAM2D(I,J)=0.0 DO 380 K=1,KBM1 DO 380 J=1,JM DO 380 I=1,IM AAM2D(I,J)=AAM2D(I,J)+AAM(I,J,K)*DZ(K) 380 CONTINUE C ELSE C DO 355 J=1,JM DO 355 I=1,IM AAM2D(I,J)=HORCON*VARHF(I,J)*TPS(I,J)/FMIN*FSM(I,J) 355 CONTINUE C ENDIF C DO 360 J=1,JM DO 360 I=1,IM KM(I,J,1)=0.0 KM(I,J,KB)=0.0 KH(I,J,1)=0.0 KH(I,J,KB)=0.0 KQ(I,J,1)=0.0 KQ(I,J,KB)=0.0 L (I,J,1)=0.0 L (I,J,KB)=0.0 Q2(I,J,1)=0.0 Q2(I,J,KB)=0.0 Q2B(I,J,1)=0.0 Q2B(I,J,KB)=0.0 Q2L(I,J,1)=0.0 Q2L(I,J,KB)=0.0 Q2LB(I,J,1)=0.0 Q2LB(I,J,KB)=0.0 360 CONTINUE C DO 370 K=2,KBM1 DO 370 J=1,JM DO 370 I=1,IM L(I,J,K) =1.*FSM(I,J) Q2B(I,J,K) =1.E-12*FSM(I,J) Q2(I,J,K) =Q2B(I,J,K) Q2LB(I,J,K)=Q2B(I,J,K)*L(I,J,K) Q2L(I,J,K) =Q2LB(I,J,K) KM(I,J,K) =UMOL*VARUF(I,J)*FSM(I,J) KQ(I,J,K) =UMOL*VARUF(I,J)*FSM(I,J) KH(I,J,K) =UMOL*VARUF(I,J)/VPRNU*FSM(I,J) 370 CONTINUE C IF(VERTMIX.EQ.'CONSTANT ') UMOL=0.0 C C********************************************************************* C C FOR PARTICLE TRACKING C IF (PARTICLE.EQ.'INCLUDE') THEN C C OPEN CORNER LOCATIONS FILE C C CORNER LOCATION CONVENTION: XCOR(I,J) = x(i-1/2, j-1/2) C (LOWER LEFT-HAND CORNER) YCOR(I,J) = y(i-1/2, j-1/2) C OPEN (UNIT=33,FILE='corner_loc',FORM='FORMATTED') C DO 1300 N=1,1000000 READ (33,*,END=1310)I,J,XCOR(I,J),YCOR(I,J) 1300 CONTINUE C C CALC. H1 AND H2 AT ELEMENT INTERFACES C 1310 DO 1320 J=1,JM DO 1320 I=1,IM H1P(I,J)=SQRT((XCOR(I+1,J)-XCOR(I,J))**2.+ + (YCOR(I+1,J)-YCOR(I,J))**2.) C H2P(I,J)=SQRT((XCOR(I,J+1)-XCOR(I,J))**2.+ + (YCOR(I,J+1)-YCOR(I,J))**2.) 1320 CONTINUE ENDIF C C********************************************************************* C C RETURN END