pkg/oasim/oasim_init_fixed.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_INIT_FIXED
                0006
                0007 C     !INTERFACE:
                0008       SUBROUTINE OASIM_INIT_FIXED(myThid)
                0009
                0010 C     !DESCRIPTION:
                0011 C     Computes constants for global irradiance calculations, reads in
                0012 C     required data files, and otherwise obtains one-time-only information
                0013 C     necessary for the run.
                0014 C
                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_INTERNAL.h"
                0023 #ifdef ALLOW_CAL
                0024 #include "cal.h"
                0025 #endif
                0026
                0027 C     !INPUT PARAMETERS:
                0028       INTEGER myThid
                0029 CEOP
                0030
                0031 #ifdef ALLOW_OASIM
                0032
                0033 C     !LOCAL VARIABLES:
                0034       INTEGER l
                0035       _RL pi1,h,c,hc,oavo,hcoavo,rlamm
                0036
                0037       _BEGIN_MASTER( myThid )
                0038
                0039 c  Degrees to radians conversion
                0040       pi1 = ACOS(-1.0D0)
                0041       pi2 = pi1*2.0 _d 0
                0042       rad = 180.0D0/pi1
                0043 c
                0044 c  Obtain Light data
                0045       CALL OASIM_INIT_ATMOW( myThid )
                0046 c
                0047 c  Quanta conversion W to uEin/s
                0048 c     Plancks constant J sec
                0049       h = 6.6256 _d -34
                0050 c     speed of light m/sec
                0051       c = 2.998 _d 8
                0052       hc = 1.0 _d 0/(h*c)
                0053 c     1/Avogadros number
                0054       oavo = 1.0 _d 0/6.023 _d 23
                0055       hcoavo = hc*oavo
                0056       DO l = 1,nlt
                0057 c      lambda in m
                0058        rlamm = FLOAT(lam(l))*1.0 _d -9
                0059 c      Watts to quanta conversion
                0060        WtouEins(l) = 1.0 _d 6*rlamm*hcoavo
                0061       ENDDO
                0062 c
                0063       o24 = 1. _d 0/24 _d 0
                0064
                0065       CALL OASIM_INIT_CLEAR( myThid )
                0066
                0067 C for aerosol data
                0068       lamaer( 1) =  250
                0069       lamaer( 2) =  325
                0070       lamaer( 3) =  350
                0071       lamaer( 4) =  375
                0072       lamaer( 5) =  400
                0073       lamaer( 6) =  425
                0074       lamaer( 7) =  450
                0075       lamaer( 8) =  475
                0076       lamaer( 9) =  500
                0077       lamaer(10) =  525
                0078       lamaer(11) =  550
                0079       lamaer(12) =  575
                0080       lamaer(13) =  600
                0081       lamaer(14) =  625
                0082       lamaer(15) =  650
                0083       lamaer(16) =  675
                0084       lamaer(17) =  700
                0085       lamaer(18) =  725
                0086       lamaer(19) =  775
                0087       lamaer(20) =  850
                0088       lamaer(21) =  950
                0089       lamaer(22) = 1050
                0090       lamaer(23) = 1150
                0091       lamaer(24) = 1250
                0092       lamaer(25) = 1350
                0093       lamaer(26) = 1450
                0094       lamaer(27) = 1550
                0095       lamaer(28) = 1650
                0096       lamaer(29) = 1750
                0097       lamaer(30) = 1900
                0098       lamaer(31) = 2200
                0099       lamaer(32) = 2900
                0100       lamaer(33) = 3700
                0101       CALL OASIM_INIT_AER( myThid )
                0102
                0103 C for light:
                0104       ozfac1 = 44.0 _d 0/6370.0 _d 0
                0105       ozfac2 = 1.0 _d 0 + 22.0 _d 0/6370.0 _d 0
                0106       p0 = 1013.25 _d 0
                0107
                0108 C for cloud model
                0109       CALL OASIM_INIT_SLINGO( myThid )
                0110
                0111       _END_MASTER( myThid )
                0112       _BARRIER
                0113
                0114       CALL OASIM_INIT_VEC( myThid )
                0115
                0116 #ifndef OASIM_READ_UNFORMATTED
                0117       CALL OASIM_EXF_INIT_FIXED( myThid )
                0118 #endif
                0119
                0120 #ifdef ALLOW_DIAGNOSTICS
                0121       IF ( useDiagnostics ) THEN
                0122         CALL OASIM_DIAGNOSTICS_INIT( myThid )
                0123       ENDIF
                0124 #endif
                0125
                0126 #endif /* ALLOW_OASIM */
                0127
                0128       RETURN
                0129       END
                0130
                0131
                0132 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0133 CBOP 0
                0134 C     !ROUTINE: OASIM_INIT_ATMOW
                0135
                0136 C     !INTERFACE:
                0137       SUBROUTINE OASIM_INIT_ATMOW( myThid)
                0138
                0139 C     !DESCRIPTION:
                0140 c     Reads in radiative transfer data: specifically atmospheric data
                0141 c     (Extraterrestrial irradiance and atmospheric optical data), and
                0142 c     water data (seawater absorption and total scattering coefficients,
                0143 c     and chl-specific absorption and total scattering data for
                0144 c     several phytoplankton groups).  PAR (350-700) begins at index 3,
                0145 c     and ends at index 17.
                0146
                0147 C     !USES:
                0148       IMPLICIT NONE
                0149 #include "SIZE.h"
                0150 #include "EEPARAMS.h"
                0151 #include "OASIM_SIZE.h"
                0152 #include "OASIM_INTERNAL.h"
                0153 #include "OASIM_PARAMS.h"
                0154 #include "OASIM_FIELDS.h"
                0155
                0156 C     !INPUT PARAMETERS:
                0157       INTEGER myThid
                0158 CEOP
                0159
                0160 #ifdef ALLOW_OASIM
                0161
                0162 C     !LOCAL VARIABLES:
                0163       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0164       character*80 title
                0165       INTEGER i, l, ilam, iUnit
                0166       _RL sFobar,sthray,soza,sawv,sao,saco2
                0167       _RL lambda,saw,sbw,sac,sbc
                0168 c
                0169 c  Atmospheric data file, atmo25b.dat
                0170       CALL MDSFINDUNIT(iUnit, myThid)
                0171       OPEN(iUnit,FILE=oasim_atmoFile,STATUS='old',FORM='formatted')
                0172       READ(iUnit,'(a80)')title
                0173       READ(iUnit,'(a80)')title
                0174       READ(iUnit,'(a80)')title
                0175       READ(iUnit,'(a80)')title
                0176       WRITE(msgBuf,'(2A)') 'OASIM_INIT_FIXED: ',
                0177      &      'ilam Fobar thray oza awv ao aco2'
                0178       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0179      &                    SQUEEZE_RIGHT , myThid)
                0180       DO l = 1,nlt
                0181        READ(iUnit,10)ilam,sFobar,sthray,soza,sawv,sao,saco2
                0182        WRITE(msgBuf,'(A,i5,6f11.4)') 'OASIM_INIT_FIXED: ',
                0183      &       ilam,sFobar,sthray,soza,sawv,sao,saco2
                0184        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0185      &                     SQUEEZE_RIGHT , myThid)
                0186        lam(l) = ilam
                0187        Fobar(l) = sFobar
                0188        thray(l) = sthray
                0189        oza(l) = soza
                0190        awv(l) = sawv
                0191        ao(l) = sao
                0192        aco2(l) = saco2
                0193        OASIM_lam(l) = ilam
                0194       ENDDO
                0195       CLOSE(iUnit)
                0196 c
                0197 10    FORMAT(i5,6f11.4)
                0198 c
                0199 c
                0200 c  Water data files, abw25b.dat
                0201       OPEN(iUnit,FILE=oasim_waterFile,STATUS='old',FORM='formatted')
                0202       DO i = 1,6
                0203        READ(iUnit,'(a50)')title
                0204 c       write(6,'(a50)')title
                0205       ENDDO
                0206       DO l = 1,nlt
                0207        READ(iUnit,20)ilam,saw,sbw
                0208 c       write(6,20)lambda,saw,sbw
                0209        IF (ilam .NE. lam(l)) THEN
                0210         WRITE(msgBuf,'(3A)') 'OASIM_INIT_FIXED:',
                0211      &   ' wavebands in OASIM_wasterFile do not match those in',
                0212      &   ' OASIM_atmoFile.'
                0213         CALL PRINT_ERROR( msgBuf , myThid)
                0214         CALL ALL_PROC_DIE( 0 )
                0215         STOP 'ABNORMAL END: S/R OASIM_INIT_FIXED'
                0216        ENDIF
                0217        aw(l) = saw
                0218        bw(l) = sbw
                0219       ENDDO
                0220       CLOSE(iUnit)
                0221 20    FORMAT(i5,f15.4,f10.4)
                0222
                0223 #endif /* ALLOW_OASIM */
                0224
                0225       RETURN
                0226       END
                0227
                0228
                0229 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0230 CBOP 0
                0231 C     !ROUTINE: OASIM_INIT_VEC
                0232
                0233 C     !INTERFACE:
                0234       SUBROUTINE OASIM_INIT_VEC( myThid )
                0235
                0236 C     !DESCRIPTION:
                0237 C  Create arrays of up, north, and east vectors in earth-fixed coords
                0238 C  for OASIM grid locations.
                0239 C
                0240 C     !USES:
                0241       IMPLICIT NONE
                0242 #include "SIZE.h"
                0243 #include "EEPARAMS.h"
                0244 #include "GRID.h"
                0245 #include "OASIM_SIZE.h"
                0246 #include "OASIM_PARAMS.h"
                0247 #include "OASIM_FIELDS.h"
                0248 #include "OASIM_INTERNAL.h"
                0249
                0250 C     !INPUT PARAMETERS:
                0251 C     myThid   :: my Thread Id number
                0252       INTEGER myThid
                0253 CEOP
                0254
                0255 #ifdef ALLOW_OASIM
                0256
                0257 C     !LOCAL VARIABLES:
                0258       INTEGER i,j,bi,bj,nv
                0259       _RL xlon,ylat,rlon,rlat,cosx,cosy,sinx,siny,upxy
                0260       _RL up(3), ea(3)
                0261
                0262 c  Convert geodetic lat/lon to Earth-centered, earth-fixed (ECEF)
                0263 c  vector (geodetic unit vector)
                0264       DO bj=myByLo(myThid),myByHi(myThid)
                0265       DO bi=myBxLo(myThid),myBxHi(myThid)
                0266       DO j=1,sNy
                0267       DO i=1,sNx
                0268         ylat = YC(i,j,bi,bj)
                0269         rlat = ylat/rad
                0270         cosy = COS(rlat)
                0271         siny = SIN(rlat)
                0272         IF (oasim_fixedLon .NE. UNSET_RL) THEN
                0273           xlon = oasim_fixedLon
                0274         ELSE
                0275           xlon = XC(i,j,bi,bj)
                0276         ENDIF
                0277         rlon = xlon/rad
                0278         cosx = COS(rlon)
                0279         sinx = SIN(rlon)
                0280         up(1) = cosy*cosx
                0281         up(2) = cosy*sinx
                0282         up(3) = siny
                0283 c
                0284 c   Compute the local East and North unit vectors
                0285         upxy = SQRT(up(1)*up(1)+up(2)*up(2))
                0286         ea(1) = -up(2)/upxy
                0287         ea(2) = up(1)/upxy
                0288         ea(3) = 0.0 _d 0
                0289 c       cross product
                0290         OASIM_no(i,j,bi,bj,1) = up(2)*ea(3) - up(3)*ea(2)
                0291         OASIM_no(i,j,bi,bj,2) = up(3)*ea(1) - up(1)*ea(3)
                0292         OASIM_no(i,j,bi,bj,3) = up(1)*ea(2) - up(2)*ea(1)
                0293         OASIM_up(i,j,bi,bj,1) = up(1)
                0294         OASIM_up(i,j,bi,bj,2) = up(2)
                0295         OASIM_up(i,j,bi,bj,3) = up(3)
                0296         OASIM_ea(i,j,bi,bj,1) = ea(1)
                0297         OASIM_ea(i,j,bi,bj,2) = ea(2)
                0298         OASIM_ea(i,j,bi,bj,3) = ea(3)
                0299       ENDDO
                0300       ENDDO
                0301       ENDDO
                0302       ENDDO
                0303
                0304 #endif /* ALLOW_OASIM */
                0305
                0306       RETURN
                0307       END
                0308
                0309
                0310 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0311 CBOP 0
                0312 C     !ROUTINE: OASIM_INIT_CLEAR
                0313
                0314 C     !INTERFACE:
                0315       SUBROUTINE OASIM_INIT_CLEAR(myThid)
                0316
                0317 C     !DESCRIPTION:
                0318
                0319 C     !USES:
                0320       IMPLICIT NONE
                0321 #include "SIZE.h"
                0322 #include "OASIM_SIZE.h"
                0323 #include "OASIM_INTERNAL.h"
                0324
                0325 C     !INPUT PARAMETERS:
                0326       INTEGER myThid
                0327 CEOP
                0328
                0329 #ifdef ALLOW_OASIM
                0330
                0331 C     !LOCAL VARIABLES:
                0332       INTEGER l
                0333       _RL fac, rlam, t, tlog
                0334       _RL a0,a1,a2,a3,b0,b1,b2,b3
                0335       DATA a0,a1 /0.9976 _d 0, 0.2194 _d 0/
                0336       DATA a2,a3 /5.554 _d -2, 6.7 _d -3/
                0337       DATA b0,b1 /5.026 _d 0, -0.01138 _d 0/
                0338       DATA b2,b3 /9.552 _d -6, -2.698 _d -9/
                0339
                0340 C for clrtrans
                0341       DO l = 1,nlt
                0342         rlamu(l) = FLOAT(lam(l))*1.0 _d -3    !lambda in um
                0343       ENDDO
                0344
                0345 C for sfcrfl
                0346       rn = 1.341 _d 0    !index of refraction of pure seawater
                0347       roair = 1.2 _d 3     !density of air g/m3
                0348       DO l = 1,nlt
                0349         rlam = FLOAT(lam(l))
                0350         IF (lam(l) .LT. 900)THEN
                0351          t = EXP(-(aw(l)+0.5 _d 0*bw(l)))
                0352          tlog = LOG(1.0 _d -36+t)
                0353          fac = a0 + a1*tlog + a2*tlog*tlog + a3*tlog*tlog*tlog
                0354          wfac(l) = MIN(fac,1.0 _d 0)
                0355          wfac(l) = MAX(fac,0.0 _d 0)
                0356         ELSE
                0357          fac = b0 + b1*rlam + b2*rlam*rlam + b3*rlam*rlam*rlam
                0358          wfac(l) = MAX(fac,0.0 _d 0)
                0359         ENDIF
                0360       ENDDO
                0361
                0362 #endif /* ALLOW_OASIM */
                0363
                0364       RETURN
                0365       END
                0366
                0367
                0368 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0369 CBOP 0
                0370 C     !ROUTINE: OASIM_INIT_AER
                0371
                0372 C     !INTERFACE:
                0373       SUBROUTINE OASIM_INIT_AER(myThid)
                0374
                0375 C     !DESCRIPTION:
                0376
                0377 C     !USES:
                0378       IMPLICIT NONE
                0379 #include "SIZE.h"
                0380 #include "OASIM_SIZE.h"
                0381 #include "OASIM_INTERNAL.h"
                0382
                0383 C     !INPUT PARAMETERS:
                0384       INTEGER myThid
                0385 CEOP
                0386
                0387 #ifdef ALLOW_OASIM
                0388
                0389 C     !LOCAL VARIABLES:
                0390       INTEGER l, laer
                0391
                0392 c     Indices to relate aerosol parameters to other light parameters
                0393       DO l = 1,nlt
                0394         DO laer = 2,nltaer
                0395          IF (lam(l) .LE. lamaer(laer))THEN
                0396           iaer(l) = laer
                0397           waer(l) = (lamaer(laer)-lam(l))/(lamaer(laer)-lamaer(laer-1))
                0398           GO TO 5
                0399          ENDIF
                0400         ENDDO
                0401 5       CONTINUE
                0402       ENDDO
                0403
                0404 #endif /* ALLOW_OASIM */
                0405
                0406       RETURN
                0407       END
                0408
                0409
                0410 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0411 CBOP 0
                0412 C     !ROUTINE: OASIM_INIT_SLINGO
                0413
                0414 C     !INTERFACE:
                0415       SUBROUTINE OASIM_INIT_SLINGO(myThid)
                0416
                0417 C     !DESCRIPTION:
                0418 C     read cloud parameters and compute
                0419 C
                0420 c       ica :: index for translating cloud arrays to light arrays
                0421
                0422 C     !USES:
                0423       IMPLICIT NONE
                0424 #include "SIZE.h"
                0425 #include "OASIM_SIZE.h"
                0426 #include "OASIM_INTERNAL.h"
                0427 #include "OASIM_SLINGO.h"
                0428
                0429 C     !INPUT PARAMETERS:
                0430       INTEGER myThid
                0431 CEOP
                0432
                0433 #ifdef ALLOW_OASIM
                0434
                0435 C     !LOCAL VARIABLES:
                0436       INTEGER l, nc, iUnit
                0437       _RL lamcld
                0438
                0439       U1 = 7.0 _d 0/4.0 _d 0
                0440       CALL MDSFINDUNIT(iUnit, myThid)
                0441       CALL OASIM_RDSLINGO(iUnit)    !reads cloud parameter file
                0442 c     Indices to relate cloud parameters to other light parameters
                0443       DO l = 1,nlt
                0444         DO nc = 1,ncld
                0445          lamcld = NINT(rnl2(nc)*1000.0 _d 0)
                0446          IF (lam(l) .LT. lamcld)THEN
                0447           ica(l) = nc
                0448           GO TO 5
                0449          ENDIF
                0450         ENDDO
                0451 5       CONTINUE
                0452       ENDDO
                0453 c       ica(nlt) = ncld
                0454
                0455 #endif /* ALLOW_OASIM */
                0456
                0457       RETURN
                0458       END
                0459
                0460
                0461 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0462 CBOP 0
                0463 C     !ROUTINE: OASIM_RDSLINGO
                0464
                0465 C     !INTERFACE:
                0466       SUBROUTINE OASIM_RDSLINGO(iUnit)
                0467
                0468 C     !DESCRIPTION:
                0469 C     Reads cloud parameters by Slingo (1989).
                0470
                0471 C     !USES:
                0472       IMPLICIT NONE
                0473 #include "SIZE.h"
                0474 #include "EEPARAMS.h"
                0475 #include "OASIM_SIZE.h"
                0476 #include "OASIM_SLINGO.h"
                0477 #include "OASIM_PARAMS.h"
                0478
                0479 C     !INPUT PARAMETERS:
                0480       INTEGER iUnit
                0481 CEOP
                0482
                0483 #ifdef ALLOW_OASIM
                0484
                0485 C     !LOCAL VARIABLES:
                0486       CHARACTER*50 title
                0487       INTEGER n
                0488       REAL*4 rn1,rn2,a4,b4,c4,d4,e4,f4
                0489 c
                0490       OPEN(iUnit,FILE=oasim_slingoFile,STATUS='old',FORM='formatted')
                0491       DO n = 1,3
                0492        READ(iUnit,'(a50)')title
                0493       ENDDO
                0494       DO n = 1,ncld
                0495        READ(iUnit,10)rn1,rn2,a4,b4,e4,f4,c4,d4
                0496        rnl1(n) = rn1
                0497        rnl2(n) = rn2
                0498        asl(n) = a4*0.01 _d 0
                0499        bsl(n) = b4
                0500        csl(n) = c4
                0501        dsl(n) = d4
                0502        esl(n) = e4
                0503        fsl(n) = f4*0.001 _d 0
                0504       ENDDO
                0505       CLOSE(iUnit)
                0506 c
                0507 10    FORMAT(2f5.2,3x,2f6.3,2f6.3,1pe9.2,1pe8.2)
                0508
                0509 #endif /* ALLOW_OASIM */
                0510
                0511       RETURN
                0512       END
                0513 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|