pkg/oasim/oasim_calc_solz.F

87dd4f7d5f Oliv*0001 #include "OASIM_OPTIONS.h"
                0002
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP 0
                0005 C     !ROUTINE: OASIM_CALC_SOLZ
                0006
                0007 C     !INTERFACE:
                0008       SUBROUTINE OASIM_CALC_SOLZ( myTime, myIter, myThid )
                0009
                0010 C     !DESCRIPTION:
                0011 C     Compute solar zenith angles
                0012 C
                0013 C     NOTE: we assume time steps do not span days, i.e. midnight is between times steps
                0014
                0015 C     !USES:
                0016       IMPLICIT NONE
                0017 #include "SIZE.h"
                0018 #include "EEPARAMS.h"
                0019 #include "PARAMS.h"
                0020 #include "OASIM_SIZE.h"
                0021 #include "OASIM_PARAMS.h"
                0022 #include "OASIM_FIELDS.h"
                0023
                0024 C     !INPUT PARAMETERS:
                0025 C     myTime   :: Current time of simulation ( s )
                0026 C     myIter   :: Current iteration number in simulation
                0027 C     myThid   :: my Thread Id number
                0028       _RL     myTime
                0029       INTEGER myIter, myThid
                0030 CEOP
                0031
                0032 #ifdef ALLOW_OASIM
                0033
                0034 C     !LOCAL VARIABLES:
                0035       INTEGER year, imon, iday, secs, lp, wd, currentdate(4)
                0036       INTEGER yday, yymmdd, difftime(4), yearstart(4)
                0037       INTEGER i, j, bi, bj
                0038       _RL t1, t2
                0039
                0040       IF (OASIM_fixedSolz .GE. 0 _d 0) THEN
                0041
                0042        DO bj=myByLo(myThid),myByHi(myThid)
                0043        DO bi=myBxLo(myThid),myBxHi(myThid)
                0044         DO j=1,sNy
                0045         DO i=1,sNx
                0046          OASIM_solz(i,j,bi,bj) = OASIM_fixedSolz
                0047         ENDDO
                0048         ENDDO
                0049        ENDDO
                0050        ENDDO
                0051
                0052       ELSE
                0053
                0054 #ifdef ALLOW_CAL
                0055        CALL cal_GetDate( myiter, mytime, currentdate, mythid )
                0056        CALL cal_convDate( currentdate,year,imon,iday,secs,lp,wd,mythid )
                0057        t1 = secs
                0058        t2 = t1 + deltaT
                0059 C      compute yearday
                0060        yymmdd = year*10000 + 101
                0061        CALL cal_fullDate( yymmdd, 0, yearstart, mythid )
                0062        CALL cal_timePassed( yearstart, currentdate, difftime, mythid )
                0063        yday = difftime(1) + 1
                0064 #else
                0065 c      make up a yearday
                0066        year = FLOOR(myTime/31104000. _d 0)
                0067        t1 = (myTime - 31104000. _d 0*year)*365. _d 0/360. _d 0
                0068        year = oasim_startYear + year
                0069 c      day of year
                0070        yday = FLOOR(t1/86400. _d 0)
                0071 c      time of day
                0072        t1 = t1 - 86400. _d 0*yday
                0073        t2 = t1 + deltaT
                0074        yday = yday + 1
                0075 #endif
                0076
                0077 #ifdef OASIM_INST_ZENITH_ANGLE
                0078        CALL OASIM_SFCSOLZ_INST(year,yday,t1,myThid)
                0079 #else
                0080 c compute (cos)average solar zenith angle from t1 to t2
                0081        CALL OASIM_SFCSOLZ(year,yday,t1,t2,oasim_dTsolz,myThid)
                0082 #endif
                0083
                0084 C     OASIM_fixedSolz
                0085       ENDIF
                0086
                0087 #endif /* ALLOW_OASIM */
                0088
                0089       RETURN
                0090       END
                0091
                0092
                0093 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0094 CBOP
                0095 C     !ROUTINE: OASIM_SFCSOLZ
                0096
                0097 C     !INTERFACE:
                0098       SUBROUTINE OASIM_SFCSOLZ(iyr,iday,t1,t2,dt,myThid)
                0099
                0100 C     !DESCRIPTION:
                0101 C     Computes solar zenith angle (stored in common block) at surface
                0102 C     given location and time.  The cosine of the zenith angle is
                0103 C     averaged over the interval [t1,t2] using the trapezoidal rule with
                0104 C     resolution dt.
                0105
                0106 C     !USES:
                0107       IMPLICIT NONE
                0108 #include "SIZE.h"
                0109 #include "EEPARAMS.h"
                0110 #include "OASIM_SIZE.h"
                0111 #include "OASIM_INTERNAL.h"
                0112 #include "OASIM_FIELDS.h"
                0113
                0114 C     !INPUT PARAMETERS:
                0115       INTEGER iyr, iday, myThid
                0116       _RL t1, t2, dt
                0117 CEOP
                0118
                0119 #ifdef ALLOW_OASIM
                0120
                0121 C     !LOCAL VARIABLES:
                0122       INTEGER i, j, bi, bj, jj, nstps, it
                0123 #ifdef OASIM_SCALAR_ZENITH_ANGLE_COMP
                0124       _RL csza
                0125 #else
                0126       _RL csza(sNx,sNy)
                0127 #endif
                0128       _RL gmt, sunz, cosunz, rsza, sza
                0129       _RL rs, sunvec(3), sunv, sunn, sune
                0130
                0131       nstps = NINT((t2-t1)/dt) + 1
                0132
                0133 c  Integrate to obtain mean cosine solar zenith angle
                0134 #ifdef OASIM_SCALAR_ZENITH_ANGLE_COMP
                0135
                0136       DO bj=myByLo(myThid),myByHi(myThid)
                0137       DO bi=myBxLo(myThid),myBxHi(myThid)
                0138       DO j=1,sNy
                0139       DO i=1,sNx
                0140        csza = 0.0 _d 0
                0141        DO it = 1,nstps,(nstps-1)
                0142         gmt = (t1 + (it-1)*dt)/3600. _d 0
                0143         CALL OASIM_SUN_VECTOR(iday,iyr,gmt,sunvec,rs)
                0144 c  Compute components of spacecraft and sun vector in the
                0145 c  vertical (up), North (no), and East (ea) vectors frame
                0146         sunv = 0.0 _d 0
                0147         sunn = 0.0 _d 0
                0148         sune = 0.0 _d 0
                0149         DO jj = 1,3
                0150          sunv = sunv + sunvec(jj)*OASIM_up(i,j,bi,bj,jj)
                0151          sunn = sunn + sunvec(jj)*OASIM_no(i,j,bi,bj,jj)
                0152          sune = sune + sunvec(jj)*OASIM_ea(i,j,bi,bj,jj)
                0153         ENDDO
                0154         sunz = rad*ATAN2(SQRT(sunn*sunn+sune*sune),sunv)
                0155         csza = csza + 0.5 _d 0*COS(sunz/rad)
                0156        ENDDO    !end of it loop
                0157        DO it = 2,nstps-1
                0158         gmt = (t1 + (it-1)*dt)/3600. _d 0
                0159         CALL OASIM_SUN_VECTOR(iday,iyr,gmt,sunvec,rs)
                0160 c  Compute components of spacecraft and sun vector in the
                0161 c  vertical (up), North (no), and East (ea) vectors frame
                0162         sunv = 0.0 _d 0
                0163         sunn = 0.0 _d 0
                0164         sune = 0.0 _d 0
                0165         DO jj = 1,3
                0166          sunv = sunv + sunvec(jj)*OASIM_up(i,j,bi,bj,jj)
                0167          sunn = sunn + sunvec(jj)*OASIM_no(i,j,bi,bj,jj)
                0168          sune = sune + sunvec(jj)*OASIM_ea(i,j,bi,bj,jj)
                0169         ENDDO
                0170         sunz = rad*ATAN2(SQRT(sunn*sunn+sune*sune),sunv)
                0171         csza = csza + COS(sunz/rad)
                0172        ENDDO   !end of it loop
                0173        cosunz = csza*rad/(nstps-1)
                0174        rsza = ACOS(cosunz/rad)
                0175        sza = rsza*rad
                0176        sza = min(sza, 90.0 _d 0)
                0177        OASIM_solz(i,j,bi,bj) = max(sza, 0.0 _d 0)
                0178       ENDDO
                0179       ENDDO
                0180       ENDDO
                0181       ENDDO
                0182
                0183 #else /* not OASIM_SCALAR_ZENITH_ANGLE_COMP */
                0184
                0185       DO bj=myByLo(myThid),myByHi(myThid)
                0186       DO bi=myBxLo(myThid),myBxHi(myThid)
                0187        DO j=1,sNy
                0188        DO i=1,sNx
                0189         csza(i,j) = 0.0 _d 0
                0190        ENDDO
                0191        ENDDO
                0192
                0193        DO it = 1,nstps,(nstps-1)
                0194         gmt = (t1 + (it-1)*dt)/3600. _d 0
                0195         CALL OASIM_SUN_VECTOR(iday,iyr,gmt,sunvec,rs)
                0196         DO j=1,sNy
                0197         DO i=1,sNx
                0198 c  Compute components of spacecraft and sun vector in the
                0199 c  vertical (up), North (no), and East (ea) vectors frame
                0200          sunv = 0.0 _d 0
                0201          sunn = 0.0 _d 0
                0202          sune = 0.0 _d 0
                0203          DO jj = 1,3
                0204           sunv = sunv + sunvec(jj)*OASIM_up(i,j,bi,bj,jj)
                0205           sunn = sunn + sunvec(jj)*OASIM_no(i,j,bi,bj,jj)
                0206           sune = sune + sunvec(jj)*OASIM_ea(i,j,bi,bj,jj)
                0207          ENDDO
                0208          sunz = rad*ATAN2(SQRT(sunn*sunn+sune*sune),sunv)
                0209          csza(i,j) = csza(i,j) + 0.5 _d 0*COS(sunz/rad)
                0210         ENDDO
                0211         ENDDO
                0212        ENDDO    !end of it loop
                0213
                0214        DO it = 2,nstps-1
                0215         gmt = (t1 + (it-1)*dt)/3600. _d 0
                0216         CALL OASIM_SUN_VECTOR(iday,iyr,gmt,sunvec,rs)
                0217         DO j=1,sNy
                0218         DO i=1,sNx
                0219 c  Compute components of spacecraft and sun vector in the
                0220 c  vertical (up), North (no), and East (ea) vectors frame
                0221          sunv = 0.0 _d 0
                0222          sunn = 0.0 _d 0
                0223          sune = 0.0 _d 0
                0224          DO jj = 1,3
                0225           sunv = sunv + sunvec(jj)*OASIM_up(i,j,bi,bj,jj)
                0226           sunn = sunn + sunvec(jj)*OASIM_no(i,j,bi,bj,jj)
                0227           sune = sune + sunvec(jj)*OASIM_ea(i,j,bi,bj,jj)
                0228          ENDDO
                0229          sunz = rad*ATAN2(SQRT(sunn*sunn+sune*sune),sunv)
                0230          csza(i,j) = csza(i,j) + COS(sunz/rad)
                0231         ENDDO
                0232         ENDDO
                0233        ENDDO   !end of it loop
                0234
                0235        DO j=1,sNy
                0236        DO i=1,sNx
                0237         cosunz = csza(i,j)/(nstps-1)
                0238         rsza = ACOS(cosunz)
                0239         sza = rsza*rad
                0240         sza = MIN(sza, 90.0 _d 0)
                0241         OASIM_solz(i,j,bi,bj) = MAX(sza, 0.0 _d 0)
                0242        ENDDO
                0243        ENDDO
                0244       ENDDO
                0245       ENDDO
                0246
                0247 #endif /* OASIM_SCALAR_ZENITH_ANGLE_COMP */
                0248
                0249 #endif /* ALLOW_OASIM */
                0250
                0251       RETURN
                0252       END
                0253
                0254 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0255 CBOP
                0256 C     !ROUTINE: OASIM_SFCSOLZ_INST
                0257
                0258 C     !INTERFACE:
                0259       SUBROUTINE OASIM_SFCSOLZ_INST(iyr,iday,t1,myThid)
                0260
                0261 C     !DESCRIPTION:
                0262 C     Computes solar zenith angle (stored in common block) at surface
                0263 C     given location and time.
                0264
                0265 C     !USES:
                0266       IMPLICIT NONE
                0267 #include "SIZE.h"
                0268 #include "EEPARAMS.h"
                0269 #include "OASIM_SIZE.h"
                0270 #include "OASIM_INTERNAL.h"
                0271 #include "OASIM_FIELDS.h"
                0272
                0273 C     !INPUT PARAMETERS:
                0274       INTEGER iyr, iday, myThid
                0275       _RL t1
                0276 CEOP
                0277
                0278 #ifdef ALLOW_OASIM
                0279
                0280 C     !LOCAL VARIABLES:
                0281       INTEGER i, j, bi, bj, jj
                0282       _RL gmt, rs, sunv, sunn, sune, sza
                0283       _RL sunvec(3)
                0284
                0285       gmt = t1/3600. _d 0
                0286       CALL OASIM_SUN_VECTOR(
                0287      I                       iday,iyr,gmt,
                0288      O                       sunvec,rs )
                0289       DO bj=myByLo(myThid),myByHi(myThid)
                0290       DO bi=myBxLo(myThid),myBxHi(myThid)
                0291        DO j=1,sNy
                0292        DO i=1,sNx
                0293 c  Compute components of spacecraft and sun vector in the
                0294 c  vertical (up), North (no), and East (ea) vectors frame
                0295         sunv = 0.0 _d 0
                0296         sunn = 0.0 _d 0
                0297         sune = 0.0 _d 0
                0298         DO jj = 1,3
                0299          sunv = sunv + sunvec(jj)*OASIM_up(i,j,bi,bj,jj)
                0300          sunn = sunn + sunvec(jj)*OASIM_no(i,j,bi,bj,jj)
                0301          sune = sune + sunvec(jj)*OASIM_ea(i,j,bi,bj,jj)
                0302         ENDDO
                0303         sza = rad*ATAN2(SQRT(sunn*sunn+sune*sune),sunv)
                0304         sza = MIN(sza, 90.0 _d 0)
                0305         OASIM_solz(i,j,bi,bj) = MAX(sza, 0.0 _d 0)
                0306        ENDDO
                0307        ENDDO
                0308       ENDDO
                0309       ENDDO
                0310
                0311 #endif /* ALLOW_OASIM */
                0312
                0313       RETURN
                0314       END