share/sint.F

   1
   2       SUBROUTINE SINT(XF,                        &
   3                    ims, ime, jms, jme, icmask ,  &
   4                    its, ite, jts, jte, nf, xstag, ystag )
   5       IMPLICIT NONE
   6       INTEGER ims, ime, jms, jme, &
   7               its, ite, jts, jte
   8
   9       LOGICAL icmask( ims:ime, jms:jme )
  10       LOGICAL xstag, ystag
  11
  12       INTEGER nf, ior
  13       REAL    one12, one24, ep
  14       PARAMETER(one12=1./12.,one24=1./24.)
  15       PARAMETER(ior=2)
  16 !
  17       REAL XF(ims:ime,jms:jme,NF)
  18 !
  19       REAL Y(ims:ime,jms:jme,-IOR:IOR),    &
  20            Z(ims:ime,jms:jme,-IOR:IOR),    &
  21            F(ims:ime,jms:jme,0:1)
  22 !
  23       INTEGER I,J,II,JJ,IIM
  24       INTEGER N2STAR, N2END, N1STAR, N1END
  25 !
  26       DATA  EP/ 1.E-10/
  27
  28       REAL W(ims:ime,jms:jme),OV(ims:ime,jms:jme),UN(ims:ime,jms:jme)
  29       REAL MXM(ims:ime,jms:jme),MN(ims:ime,jms:jme)
  30       REAL FL(ims:ime,jms:jme,0:1)
  31       REAL XIG(NF*NF), XJG(NF*NF)  ! NF is parent to child grid refinement ratio
  32       integer rr
  33
  34       REAL rioff, rjoff
  35 !
  36       REAL donor, y1, y2, a
  37       DONOR(Y1,Y2,A)=(Y1*AMAX1(0.,SIGN(1.,A))-Y2*AMIN1(0.,SIGN(1.,A)))*A
  38       REAL tr4, ym1, y0, yp1, yp2
  39       TR4(YM1,Y0,YP1,YP2,A)=A*ONE12*(7.*(YP1+Y0)-(YP2+YM1))               &
  40        -A*A*ONE24*(15.*(YP1-Y0)-(YP2-YM1))-A*A*A*ONE12*((YP1+Y0)          &
  41        -(YP2+YM1))+A*A*A*A*ONE24*(3.*(YP1-Y0)-(YP2-YM1))
  42       REAL pp, pn, x
  43       PP(X)=AMAX1(0.,X)
  44       PN(X)=AMIN1(0.,X)
  45
  46       rr = nint(sqrt(float(nf)))
  47 !!      write(6,*) ' nf, rr are ',nf,rr
  48
  49       rioff = 0
  50       rjoff = 0
  51       if(xstag .and. (mod(rr,2) .eq. 0)) rioff = 1.
  52       if(ystag .and. (mod(rr,2) .eq. 0)) rjoff = 1.
  53
  54       DO I=1,rr
  55         DO J=1,rr
  56           XIG(J+(I-1)*rr)=(float(rr)-1.-rioff)/float(2*rr)-FLOAT(J-1)*1./float(rr)
  57           XJG(J+(I-1)*rr)=(float(rr)-1.-rjoff)/float(2*rr)-FLOAT(I-1)*1./float(rr)
  58         ENDDO
  59       ENDDO
  60
  61       N2STAR = jts
  62       N2END  = jte
  63       N1STAR = its
  64       N1END  = ite
  65
  66       DO 2000 IIM=1,NF
  67 !
  68 !  HERE STARTS RESIDUAL ADVECTION
  69 !
  70         DO 9000 JJ=N2STAR,N2END
  71           DO 50 J=-IOR,IOR
  72
  73             DO 51 I=-IOR,IOR
  74               DO 511 II=N1STAR,N1END
  75                 IF ( icmask(II,JJ) ) Y(II,JJ,I)=XF(II+I,JJ+J,IIM)
  76   511         CONTINUE
  77    51       CONTINUE
  78
  79             DO 811 II=N1STAR,N1END
  80               IF ( icmask(II,JJ) ) THEN
  81                 FL(II,JJ,0)=DONOR(Y(II,JJ,-1),Y(II,JJ,0),XIG(IIM))
  82                 FL(II,JJ,1)=DONOR(Y(II,JJ,0),Y(II,JJ,1),XIG(IIM))
  83               ENDIF
  84   811         CONTINUE
  85             DO 812 II=N1STAR,N1END
  86               IF ( icmask(II,JJ) ) W(II,JJ)=Y(II,JJ,0)-(FL(II,JJ,1)-FL(II,JJ,0))
  87   812         CONTINUE
  88             DO 813 II=N1STAR,N1END
  89               IF ( icmask(II,JJ) ) THEN
  90                 MXM(II,JJ)=                                             &
  91                          AMAX1(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1),       &
  92                          W(II,JJ))
  93                 MN(II,JJ)=AMIN1(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1),W(II,JJ))
  94               ENDIF
  95   813         CONTINUE
  96             DO 312 II=N1STAR,N1END
  97               IF ( icmask(II,JJ) ) THEN
  98                 F(II,JJ,0)=                                               &
  99                            TR4(Y(II,JJ,-2),Y(II,JJ,-1),Y(II,JJ,0),        &
 100                            Y(II,JJ,1),XIG(IIM))
 101                 F(II,JJ,1)=                                                 &
 102                          TR4(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1),Y(II,JJ,2),&
 103                          XIG(IIM))
 104                 ENDIF
 105   312         CONTINUE
 106             DO 822 II=N1STAR,N1END
 107               IF ( icmask(II,JJ) ) THEN
 108                 F(II,JJ,0)=F(II,JJ,0)-FL(II,JJ,0)
 109                 F(II,JJ,1)=F(II,JJ,1)-FL(II,JJ,1)
 110               ENDIF
 111   822         CONTINUE
 112             DO 823 II=N1STAR,N1END
 113               IF ( icmask(II,JJ) ) THEN
 114                 OV(II,JJ)=(MXM(II,JJ)-W(II,JJ))/(-PN(F(II,JJ,1))+         &
 115                         PP(F(II,JJ,0))+EP)
 116                 UN(II,JJ)=(W(II,JJ)-MN(II,JJ))/(PP(F(II,JJ,1))-             &
 117                       PN(F(II,JJ,0))+EP)
 118               ENDIF
 119   823         CONTINUE
 120             DO 824 II=N1STAR,N1END
 121               IF ( icmask(II,JJ) ) THEN
 122                 F(II,JJ,0)=PP(F(II,JJ,0))*AMIN1(1.,OV(II,JJ))+            &
 123                            PN(F(II,JJ,0))*AMIN1(1.,UN(II,JJ))
 124                 F(II,JJ,1)=PP(F(II,JJ,1))*AMIN1(1.,UN(II,JJ))+            &
 125                            PN(F(II,JJ,1))*AMIN1(1.,OV(II,JJ))
 126               ENDIF
 127   824         CONTINUE
 128             DO 825 II=N1STAR,N1END
 129               IF ( icmask(II,JJ) ) THEN
 130                 Y(II,JJ,0)=W(II,JJ)-(F(II,JJ,1)-F(II,JJ,0))
 131               ENDIF
 132   825         CONTINUE
 133             DO 361 II=N1STAR,N1END
 134               IF ( icmask(II,JJ) ) Z(II,JJ,J)=Y(II,JJ,0)
 135   361         CONTINUE
 136 !
 137 !  END IF FIRST J LOOP
 138 !
 139  8000       CONTINUE
 140    50     CONTINUE
 141
 142           DO 911 II=N1STAR,N1END
 143             IF ( icmask(II,JJ) ) THEN
 144               FL(II,JJ,0)=DONOR(Z(II,JJ,-1),Z(II,JJ,0),XJG(IIM))
 145               FL(II,JJ,1)=DONOR(Z(II,JJ,0),Z(II,JJ,1),XJG(IIM))
 146             ENDIF
 147   911       CONTINUE
 148           DO 912 II=N1STAR,N1END
 149             IF ( icmask(II,JJ) ) W(II,JJ)=Z(II,JJ,0)-(FL(II,JJ,1)-FL(II,JJ,0))
 150   912       CONTINUE
 151           DO 913 II=N1STAR,N1END
 152             IF ( icmask(II,JJ) ) THEN
 153               MXM(II,JJ)=AMAX1(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),W(II,JJ))
 154               MN(II,JJ)=AMIN1(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),W(II,JJ))
 155             ENDIF
 156   913       CONTINUE
 157           DO 412 II=N1STAR,N1END
 158             IF ( icmask(II,JJ) ) THEN
 159               F(II,JJ,0)=                                                 &
 160                          TR4(Z(II,JJ,-2),Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1)&
 161                          ,XJG(IIM))
 162               F(II,JJ,1)=                                                   &
 163                          TR4(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),Z(II,JJ,2),  &
 164                          XJG(IIM))
 165             ENDIF
 166   412       CONTINUE
 167           DO 922 II=N1STAR,N1END
 168             IF ( icmask(II,JJ) ) THEN
 169               F(II,JJ,0)=F(II,JJ,0)-FL(II,JJ,0)
 170               F(II,JJ,1)=F(II,JJ,1)-FL(II,JJ,1)
 171             ENDIF
 172   922       CONTINUE
 173           DO 923 II=N1STAR,N1END
 174             IF ( icmask(II,JJ) ) THEN
 175               OV(II,JJ)=(MXM(II,JJ)-W(II,JJ))/(-PN(F(II,JJ,1))+           &
 176                         PP(F(II,JJ,0))+EP)
 177               UN(II,JJ)=(W(II,JJ)-MN(II,JJ))/(PP(F(II,JJ,1))-PN(F(II,JJ,0))+ &
 178                       EP)
 179             ENDIF
 180   923       CONTINUE
 181           DO 924 II=N1STAR,N1END
 182             IF ( icmask(II,JJ) ) THEN
 183               F(II,JJ,0)=PP(F(II,JJ,0))*AMIN1(1.,OV(II,JJ))+PN(F(II,JJ,0))  &
 184                          *AMIN1(1.,UN(II,JJ))
 185               F(II,JJ,1)=PP(F(II,JJ,1))*AMIN1(1.,UN(II,JJ))+PN(F(II,JJ,1))  &
 186                          *AMIN1(1.,OV(II,JJ))
 187             ENDIF
 188   924     CONTINUE
 189  9000   CONTINUE
 190         DO 925 JJ=N2STAR,N2END
 191           DO 925 II=N1STAR,N1END
 192             IF ( icmask(II,JJ) ) XF(II,JJ,IIM)=W(II,JJ)-(F(II,JJ,1)-F(II,JJ,0))
 193   925     CONTINUE
 194
 195 !
 196  2000 CONTINUE
 197       RETURN
 198       END
 199
 200 ! Version of sint that replaces mask with detailed ranges for avoiding boundaries
 201 ! may help performance by getting the conditionals out of innner loops
 202
 203       SUBROUTINE SINTB(XF1, XF ,                  &
 204                    ims, ime, jms, jme, icmask ,  &
 205                    its, ite, jts, jte, nf, xstag, ystag )
 206       IMPLICIT NONE
 207       INTEGER ims, ime, jms, jme, &
 208               its, ite, jts, jte
 209
 210       LOGICAL icmask( ims:ime, jms:jme )
 211       LOGICAL xstag, ystag
 212
 213       INTEGER nf, ior
 214       REAL    one12, one24, ep
 215       PARAMETER(one12=1./12.,one24=1./24.)
 216       PARAMETER(ior=2)
 217 !
 218       REAL XF(ims:ime,jms:jme,NF)
 219       REAL XF1(ims:ime,jms:jme,NF)
 220 !
 221       REAL Y(-IOR:IOR),    &
 222            Z(ims:ime,-IOR:IOR),    &
 223            F(0:1)
 224 !
 225       INTEGER I,J,II,JJ,IIM
 226       INTEGER N2STAR, N2END, N1STAR, N1END
 227 !
 228       DATA  EP/ 1.E-10/
 229 !
 230 !      PARAMETER(NONOS=1)
 231 !      PARAMETER(N1OS=N1*NONOS+1-NONOS,N2OS=N2*NONOS+1-NONOS)
 232 !
 233       REAL W,OV,UN
 234       REAL MXM,MN
 235       REAL FL(0:1)
 236       REAL XIG(NF*NF), XJG(NF*NF)  ! NF is the parent to child grid refinement ratio
 237       integer rr
 238
 239       REAL rioff, rjoff
 240 !
 241       REAL donor, y1, y2, a
 242       DONOR(Y1,Y2,A)=(Y1*AMAX1(0.,SIGN(1.,A))-Y2*AMIN1(0.,SIGN(1.,A)))*A
 243       REAL tr4, ym1, y0, yp1, yp2
 244       TR4(YM1,Y0,YP1,YP2,A)=A*ONE12*(7.*(YP1+Y0)-(YP2+YM1))               &
 245        -A*A*ONE24*(15.*(YP1-Y0)-(YP2-YM1))-A*A*A*ONE12*((YP1+Y0)          &
 246        -(YP2+YM1))+A*A*A*A*ONE24*(3.*(YP1-Y0)-(YP2-YM1))
 247       REAL pp, pn, x
 248       PP(X)=AMAX1(0.,X)
 249       PN(X)=AMIN1(0.,X)
 250
 251       rr = nint(sqrt(float(nf)))
 252
 253       rioff = 0
 254       rjoff = 0
 255       if(xstag .and. (mod(rr,2) .eq. 0)) rioff = 1.
 256       if(ystag .and. (mod(rr,2) .eq. 0)) rjoff = 1.
 257
 258       DO I=1,rr
 259         DO J=1,rr
 260           XIG(J+(I-1)*rr)=(float(rr)-1.-rioff)/float(2*rr)-FLOAT(J-1)*1./float(rr)
 261           XJG(J+(I-1)*rr)=(float(rr)-1.-rjoff)/float(2*rr)-FLOAT(I-1)*1./float(rr)
 262         ENDDO
 263       ENDDO
 264
 265       N2STAR = jts
 266       N2END  = jte
 267       N1STAR = its
 268       N1END  = ite
 269
 270       DO 2000 IIM=1,NF
 271 !
 272 !  HERE STARTS RESIDUAL ADVECTION
 273 !
 274         DO 9000 JJ=N2STAR,N2END
 275 !cdir unroll=5
 276           DO 50 J=-IOR,IOR
 277
 278 !cdir unroll=5
 279               DO 511 II=N1STAR,N1END
 280                 Y(-2)=XF1(II-2,JJ+J,IIM)
 281                 Y(-1)=XF1(II-1,JJ+J,IIM)
 282                 Y(0)=XF1(II,JJ+J,IIM)
 283                 Y(1)=XF1(II+1,JJ+J,IIM)
 284                 Y(2)=XF1(II+2,JJ+J,IIM)
 285
 286               FL(0)=DONOR(Y(-1),Y(0),XIG(IIM))
 287               FL(1)=DONOR(Y(0),Y(1),XIG(IIM))
 288               W=Y(0)-(FL(1)-FL(0))
 289               MXM=                                             &
 290                        AMAX1(Y(-1),Y(0),Y(1),       &
 291                        W)
 292               MN=AMIN1(Y(-1),Y(0),Y(1),W)
 293               F(0)=                                               &
 294                    TR4(Y(-2),Y(-1),Y(0),        &
 295                    Y(1),XIG(IIM))
 296               F(1)=                                                 &
 297                        TR4(Y(-1),Y(0),Y(1),Y(2),&
 298                        XIG(IIM))
 299               F(0)=F(0)-FL(0)
 300               F(1)=F(1)-FL(1)
 301               OV=(MXM-W)/(-PN(F(1))+         &
 302                       PP(F(0))+EP)
 303               UN=(W-MN)/(PP(F(1))-             &
 304                     PN(F(0))+EP)
 305               F(0)=PP(F(0))*AMIN1(1.,OV)+            &
 306                    PN(F(0))*AMIN1(1.,UN)
 307               F(1)=PP(F(1))*AMIN1(1.,UN)+            &
 308                    PN(F(1))*AMIN1(1.,OV)
 309               Y(0)=W-(F(1)-F(0))
 310               Z(II,J)=Y(0)
 311   511         CONTINUE
 312 !
 313 !  END IF FIRST J LOOP
 314 !
 315  8000       CONTINUE
 316    50     CONTINUE
 317
 318           DO 911 II=N1STAR,N1END
 319             FL(0)=DONOR(Z(II,-1),Z(II,0),XJG(IIM))
 320             FL(1)=DONOR(Z(II,0),Z(II,1),XJG(IIM))
 321             W=Z(II,0)-(FL(1)-FL(0))
 322             MXM=AMAX1(Z(II,-1),Z(II,0),Z(II,1),W)
 323              MN=AMIN1(Z(II,-1),Z(II,0),Z(II,1),W)
 324             F(0)=                                                 &
 325                  TR4(Z(II,-2),Z(II,-1),Z(II,0),Z(II,1)&
 326                  ,XJG(IIM))
 327             F(1)=                                                   &
 328                  TR4(Z(II,-1),Z(II,0),Z(II,1),Z(II,2),  &
 329                  XJG(IIM))
 330             F(0)=F(0)-FL(0)
 331             F(1)=F(1)-FL(1)
 332             OV=(MXM-W)/(-PN(F(1))+           &
 333                PP(F(0))+EP)
 334             UN=(W-MN)/(PP(F(1))-PN(F(0))+ &
 335                     EP)
 336             F(0)=PP(F(0))*AMIN1(1.,OV)+PN(F(0))  &
 337                  *AMIN1(1.,UN)
 338             F(1)=PP(F(1))*AMIN1(1.,UN)+PN(F(1))  &
 339                        *AMIN1(1.,OV)
 340             XF(II,JJ,IIM)=W-(F(1)-F(0))
 341   911     CONTINUE
 342  9000   CONTINUE
 343
 344 !
 345  2000 CONTINUE
 346       RETURN
 347       END
 348
 349