diff --git a/MSFPS_8f.html b/MSFPS_8f.html index 973e12885..b34dd4130 100644 --- a/MSFPS_8f.html +++ b/MSFPS_8f.html @@ -3,7 +3,7 @@ - + UPP: MSFPS.f File Reference @@ -30,7 +30,7 @@
UPP -  V11.0.0 +  11.0.0
@@ -38,7 +38,7 @@ - + @@ -109,19 +109,19 @@
-

msfps() computes the map scale factor for a polar stereographic grid at a give latitude. +

msfps() computes the map scale factor for a polar stereographic grid at a give latitude. More...

Go to the source code of this file.

- - + + +

Functions/Subroutines

-subroutine MSFPS (LAT, TRUELAT1, MSF)
 
subroutine msfps (LAT, TRUELAT1, MSF)
 msfps() computes the map scale factor for a polar stereographic grid at a give latitude. More...
 

Detailed Description

-

msfps() computes the map scale factor for a polar stereographic grid at a give latitude.

+

msfps() computes the map scale factor for a polar stereographic grid at a give latitude.

This subroutine computes the map scale factor for a polar stereographic grid at a give latitude.

Parameters
@@ -142,7 +142,55 @@

Program History Log

Date
2006-11-01

Definition in file MSFPS.f.

- +

Function/Subroutine Documentation

+ +
+
+
+ + + + + + + + + + + + + + + + + + + + + + + +
subroutine msfps (real, intent(in) LAT,
real, intent(in) TRUELAT1,
real, intent(out) MSF 
)
+
+ +

msfps() computes the map scale factor for a polar stereographic grid at a give latitude.

+

This subroutine computes the map scale factor for a polar stereographic grid at a give latitude.

+
Parameters
+ + + + +
[in]LATLatitude at which map factor is valid.
[in]TRUELAT1TRUELAT 1.
[out]MSFMap scale factor.
+
+
+ +

Definition at line 26 of file MSFPS.f.

+ +

Referenced by initpost(), initpost_gfs_nems_mpiio(), and initpost_netcdf().

+ +
+
+
diff --git a/MSFPS_8f.js b/MSFPS_8f.js index a399f8b27..be5fc74b4 100644 --- a/MSFPS_8f.js +++ b/MSFPS_8f.js @@ -1,4 +1,4 @@ var MSFPS_8f = [ - [ "MSFPS", "MSFPS_8f.html#a61f85c739a82a005bcba00dd4b7443ee", null ] + [ "msfps", "MSFPS_8f.html#a6c8f1d60cbab1f64dc65f900dc623b2e", null ] ]; \ No newline at end of file diff --git a/MSFPS_8f_source.html b/MSFPS_8f_source.html index 64baf45fc..16a7a66bc 100644 --- a/MSFPS_8f_source.html +++ b/MSFPS_8f_source.html @@ -3,7 +3,7 @@ - + UPP: MSFPS.f Source File @@ -30,7 +30,7 @@
UPP -  V11.0.0 +  11.0.0
@@ -38,7 +38,7 @@ - + @@ -107,39 +107,43 @@
Go to the documentation of this file.
1 
-
16  SUBROUTINE msfps(LAT,TRUELAT1,MSF)
-
17 
-
18 
-
19 ! Computes the map scale factor for a Polar Stereographic grid at a given
-
20 ! latitude.
-
21 
-
22  IMPLICIT NONE
-
23 
-
24 ! Define some private constants
-
25 !
-
26  REAL, PARAMETER :: pi = 3.1415927
-
27  REAL, PARAMETER :: rad_per_deg = pi / 180.
+
16 
+
17 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
25 
+
26  SUBROUTINE msfps(LAT,TRUELAT1,MSF)
+
27 
28 
-
29  REAL, INTENT(IN) :: lat ! latitude where msf is requested
-
30  REAL, INTENT(IN) :: truelat1
-
31  REAL, INTENT(OUT) :: msf
-
32 
-
33  REAL :: psi1, psix, pole
-
34 
-
35  IF (truelat1 >= 0.) THEN
-
36  psi1 = (90. - truelat1) * rad_per_deg
-
37  pole =90.
-
38  ELSE
-
39  psi1 = (90. + truelat1) * rad_per_deg
-
40  pole = -90.
-
41  ENDIF
+
29 ! Computes the map scale factor for a Polar Stereographic grid at a given
+
30 ! latitude.
+
31 
+
32  IMPLICIT NONE
+
33 
+
34 ! Define some private constants
+
35 !
+
36  REAL, PARAMETER :: pi = 3.1415927
+
37  REAL, PARAMETER :: rad_per_deg = pi / 180.
+
38 
+
39  REAL, INTENT(IN) :: lat ! latitude where msf is requested
+
40  REAL, INTENT(IN) :: truelat1
+
41  REAL, INTENT(OUT) :: msf
42 
-
43  psix = (pole - lat)*rad_per_deg
-
44  msf = ((1.+cos(psi1))/(1.0 + cos(psix)))
-
45  RETURN
-
46 
-
47  END SUBROUTINE msfps
-
48 
+
43  REAL :: psi1, psix, pole
+
44 
+
45  IF (truelat1 >= 0.) THEN
+
46  psi1 = (90. - truelat1) * rad_per_deg
+
47  pole =90.
+
48  ELSE
+
49  psi1 = (90. + truelat1) * rad_per_deg
+
50  pole = -90.
+
51  ENDIF
+
52 
+
53  psix = (pole - lat)*rad_per_deg
+
54  msf = ((1.+cos(psi1))/(1.0 + cos(psix)))
+
55  RETURN
+
56 
+
57  END SUBROUTINE msfps
+
58 
+
msfps
subroutine msfps(LAT, TRUELAT1, MSF)
msfps() computes the map scale factor for a polar stereographic grid at a give latitude.
Definition: MSFPS.f:26
@@ -148,7 +152,7 @@
  • sorc
  • ncep_post.fd
  • MSFPS.f
  • Generated by - doxygen 1.8.6
    • + doxygen 1.8.5 diff --git a/NGMFLD_8f.html b/NGMFLD_8f.html index adb57ccd8..5d2249bef 100644 --- a/NGMFLD_8f.html +++ b/NGMFLD_8f.html @@ -3,7 +3,7 @@ - + UPP: NGMFLD.f File Reference @@ -30,7 +30,7 @@
    UPP -  V11.0.0 +  11.0.0
    @@ -38,7 +38,7 @@ - + @@ -109,21 +109,21 @@
    -

    ngmfld() computes layer mean NGM fields +

    ngmfld() computes layer mean NGM fields More...

    Go to the source code of this file.

    - - + + +

    Functions/Subroutines

    -subroutine NGMFLD (RH4710, RH4796, RH1847, RH8498, QM8510)
     
    subroutine ngmfld (RH4710, RH4796, RH1847, RH8498, QM8510)
     ngmfld() computes layer mean NGM fields More...
     

    Detailed Description

    -

    ngmfld() computes layer mean NGM fields

    +

    ngmfld() computes layer mean NGM fields

    This routine computes a handful of NGM layer mean fields. This is done to provide a fully complete ETA NGM look-alike output file.

    -

    The sigma (layer) fields computed bu this routine are tabulated below.

    +

    The sigma (layer) fields computed by this routine are tabulated below.

    @@ -179,7 +179,69 @@

    Program History Log

    Date
    1992-12-22

    Definition in file NGMFLD.f.

    - +

    Function/Subroutine Documentation

    + +
    +
    +
    Sigma (layer) Field(s)
    + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
    subroutine ngmfld (real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(out) RH4710,
    real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(out) RH4796,
    real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(out) RH1847,
    real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(out) RH8498,
    real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(out) QM8510 
    )
    +
    + +

    ngmfld() computes layer mean NGM fields

    +

    This routine computes a handful of NGM layer mean fields.

    +
    Parameters
    + + + + + + +
    [out]RH4710Sigma layer 0.47-1.00 mean relative humidity.
    [out]RH4796Sigma layer 0.47-0.96 mean relative humidity.
    [out]RH1847Sigma layer 0.18-0.47 mean relative humidity.
    [out]RH8498Sigma layer 0.84-0.98 mean relative humidity.
    [out]QM8510Sigma layer 0.85-1.00 mean moisture convergence.
    +
    +
    + +

    Definition at line 59 of file NGMFLD.f.

    + +

    Referenced by miscln().

    + +
    +
    +
    diff --git a/NGMFLD_8f.js b/NGMFLD_8f.js index 9ce41f3d2..eb5fca342 100644 --- a/NGMFLD_8f.js +++ b/NGMFLD_8f.js @@ -1,4 +1,4 @@ var NGMFLD_8f = [ - [ "NGMFLD", "NGMFLD_8f.html#abbd3cc6c30e2fb5149a0658ab723b97b", null ] + [ "ngmfld", "NGMFLD_8f.html#ac7aab347a5426467a81ff0fbc0e80457", null ] ]; \ No newline at end of file diff --git a/NGMFLD_8f_source.html b/NGMFLD_8f_source.html index 72079dedc..43eae01b1 100644 --- a/NGMFLD_8f_source.html +++ b/NGMFLD_8f_source.html @@ -3,7 +3,7 @@ - + UPP: NGMFLD.f Source File @@ -30,7 +30,7 @@
    UPP -  V11.0.0 +  11.0.0
    @@ -38,7 +38,7 @@ - + @@ -107,195 +107,200 @@
    Go to the documentation of this file.
    1 
    -
    46  SUBROUTINE ngmfld(RH4710,RH4796,RH1847,RH8498,QM8510)
    -
    47 
    -
    48 !
    -
    49 !
    -
    50 ! INCLUDE PARAMETERS
    -
    51  use vrbls3d, only: q, uh, vh, pint, alpint, zint, t
    -
    52  use masks, only: lmh
    -
    53  use params_mod, only: d00, d50, h1m12, pq0, a2, a3, a4, h1, d01, small
    -
    54  use ctlblk_mod, only: jsta, jend, lm, jsta_2l, jend_2u, jsta_m2, jend_m2,&
    -
    55  spval, im, &
    -
    56  ista, iend, ista_2l, iend_2u, ista_m2, iend_m2, ista_m, iend_m
    -
    57 !
    -
    58 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    -
    59  implicit none
    -
    60 !
    -
    61  real,PARAMETER :: sig100=1.00000, sig98=0.98230, sig96=0.96470
    -
    62  real,PARAMETER :: sig89 =0.89671, sig85=0.85000, sig84=0.84368
    -
    63  real,PARAMETER :: sig78 =0.78483, sig47=0.47191, sig18=0.18018
    -
    64 !
    -
    65 ! DECLARE VARIABLES.
    -
    66  LOGICAL got8510,got4710,got4796,got1847,got8498
    -
    67  REAL,dimension(ista_2l:iend_2u,jsta_2l:jend_2u),intent(out) :: qm8510,rh4710,rh8498, &
    -
    68  rh4796,rh1847
    -
    69  REAL,dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: z8510,z4710,z8498,z4796,z1847
    -
    70  real,dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: q1d, u1d, v1d, qcnvg
    -
    71 !
    -
    72  integer i,j,l
    -
    73  real p100,p85,p98,p96,p84,p47,p18,alpm,de,pm,tm,qm, &
    -
    74  qmcvg,qs,rh,dz
    -
    75 !********************************************************************
    -
    76 ! START NGMFLD HERE.
    +
    46 
    +
    47 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    +
    58 
    +
    59  SUBROUTINE ngmfld(RH4710,RH4796,RH1847,RH8498,QM8510)
    +
    60 
    +
    61 !
    +
    62 !
    +
    63 ! INCLUDE PARAMETERS
    +
    64  use vrbls3d, only: q, uh, vh, pint, alpint, zint, t
    +
    65  use masks, only: lmh
    +
    66  use params_mod, only: d00, d50, h1m12, pq0, a2, a3, a4, h1, d01, small
    +
    67  use ctlblk_mod, only: jsta, jend, lm, jsta_2l, jend_2u, jsta_m2, jend_m2,&
    +
    68  spval, im, &
    +
    69  ista, iend, ista_2l, iend_2u, ista_m2, iend_m2, ista_m, iend_m
    +
    70 !
    +
    71 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    +
    72  implicit none
    +
    73 !
    +
    74  real,PARAMETER :: sig100=1.00000, sig98=0.98230, sig96=0.96470
    +
    75  real,PARAMETER :: sig89 =0.89671, sig85=0.85000, sig84=0.84368
    +
    76  real,PARAMETER :: sig78 =0.78483, sig47=0.47191, sig18=0.18018
    77 !
    -
    78 ! INITIALIZE ARRAYS.
    -
    79 !$omp parallel do private(i,j)
    -
    80  DO j=jsta,jend
    -
    81  DO i=ista,iend
    -
    82  qm8510(i,j) = d00
    -
    83  rh4710(i,j) = d00
    -
    84  rh8498(i,j) = d00
    -
    85  rh4796(i,j) = d00
    -
    86  rh1847(i,j) = d00
    -
    87  z8510(i,j) = d00
    -
    88  z8498(i,j) = d00
    -
    89  z4710(i,j) = d00
    -
    90  z4796(i,j) = d00
    -
    91  z1847(i,j) = d00
    -
    92  ENDDO
    -
    93  ENDDO
    -
    94 !
    -
    95 ! LOOP OVER HORIZONTAL GRID.
    -
    96 !
    -
    97 !!$omp parallel do &
    -
    98 ! & private(dz,p100,p18,p47,p84,p85, &
    -
    99 ! & p96,p98,pm,qdiv,qk,qkhn,qkhs,qkm1,qm,qm8510, &
    -
    100 ! & qmcvg,qs,qudx,qvdy,r2dx,r2dy,rh,rh1847,rh4710, &
    -
    101 ! & rh4796,rh8498,tm,tmt0,tmt15,z1847,z4710,z4796, &
    -
    102 ! & z8498,z8510,q1d,u1d,v1d,qcnvg)
    -
    103 
    -
    104  DO l=1,lm
    -
    105 ! COMPUTE MOISTURE CONVERGENCE
    -
    106 !$omp parallel do private(i,j)
    -
    107  DO j=jsta_2l,jend_2u
    -
    108  DO i=ista_2l,iend_2u
    -
    109  q1d(i,j) = q(i,j,l)
    -
    110  u1d(i,j) = uh(i,j,l)
    -
    111  v1d(i,j) = vh(i,j,l)
    -
    112  ENDDO
    -
    113  ENDDO
    -
    114  CALL calmcvg(q1d,u1d,v1d,qcnvg)
    -
    115 ! COMPUTE MOISTURE CONVERGENCE
    -
    116  DO j=jsta_m2,jend_m2
    -
    117  DO i=ista_m,iend_m
    -
    118 !
    -
    119 ! SET TARGET PRESSURES.
    -
    120 
    -
    121  p100 = pint(i,j,nint(lmh(i,j)))
    -
    122  p98 = sig98*p100
    -
    123  p96 = sig96*p100
    -
    124  p85 = sig85*p100
    -
    125  p84 = sig84*p100
    -
    126  p47 = sig47*p100
    -
    127  p18 = sig18*p100
    -
    128 !
    -
    129 !
    -
    130 ! COMPUTE LAYER MEAN FIELDS AT THE GIVEN K.
    +
    78 ! DECLARE VARIABLES.
    +
    79  LOGICAL got8510,got4710,got4796,got1847,got8498
    +
    80  REAL,dimension(ista_2l:iend_2u,jsta_2l:jend_2u),intent(out) :: qm8510,rh4710,rh8498, &
    +
    81  rh4796,rh1847
    +
    82  REAL,dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: z8510,z4710,z8498,z4796,z1847
    +
    83  real,dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: q1d, u1d, v1d, qcnvg
    +
    84 !
    +
    85  integer i,j,l
    +
    86  real p100,p85,p98,p96,p84,p47,p18,alpm,de,pm,tm,qm, &
    +
    87  qmcvg,qs,rh,dz
    +
    88 !********************************************************************
    +
    89 ! START NGMFLD HERE.
    +
    90 !
    +
    91 ! INITIALIZE ARRAYS.
    +
    92 !$omp parallel do private(i,j)
    +
    93  DO j=jsta,jend
    +
    94  DO i=ista,iend
    +
    95  qm8510(i,j) = d00
    +
    96  rh4710(i,j) = d00
    +
    97  rh8498(i,j) = d00
    +
    98  rh4796(i,j) = d00
    +
    99  rh1847(i,j) = d00
    +
    100  z8510(i,j) = d00
    +
    101  z8498(i,j) = d00
    +
    102  z4710(i,j) = d00
    +
    103  z4796(i,j) = d00
    +
    104  z1847(i,j) = d00
    +
    105  ENDDO
    +
    106  ENDDO
    +
    107 !
    +
    108 ! LOOP OVER HORIZONTAL GRID.
    +
    109 !
    +
    110 !!$omp parallel do &
    +
    111 ! & private(dz,p100,p18,p47,p84,p85, &
    +
    112 ! & p96,p98,pm,qdiv,qk,qkhn,qkhs,qkm1,qm,qm8510, &
    +
    113 ! & qmcvg,qs,qudx,qvdy,r2dx,r2dy,rh,rh1847,rh4710, &
    +
    114 ! & rh4796,rh8498,tm,tmt0,tmt15,z1847,z4710,z4796, &
    +
    115 ! & z8498,z8510,q1d,u1d,v1d,qcnvg)
    +
    116 
    +
    117  DO l=1,lm
    +
    118 ! COMPUTE MOISTURE CONVERGENCE
    +
    119 !$omp parallel do private(i,j)
    +
    120  DO j=jsta_2l,jend_2u
    +
    121  DO i=ista_2l,iend_2u
    +
    122  q1d(i,j) = q(i,j,l)
    +
    123  u1d(i,j) = uh(i,j,l)
    +
    124  v1d(i,j) = vh(i,j,l)
    +
    125  ENDDO
    +
    126  ENDDO
    +
    127  CALL calmcvg(q1d,u1d,v1d,qcnvg)
    +
    128 ! COMPUTE MOISTURE CONVERGENCE
    +
    129  DO j=jsta_m2,jend_m2
    +
    130  DO i=ista_m,iend_m
    131 !
    -
    132 ! COMPUTE P, Z, T, AND Q AT THE MIDPOINT OF THE CURRENT ETA LAYER.
    -
    133  alpm = d50*(alpint(i,j,l)+alpint(i,j,l+1))
    -
    134  dz = zint(i,j,l)-zint(i,j,l+1)
    -
    135  pm = exp(alpm)
    -
    136  tm = t(i,j,l)
    -
    137  qm = q(i,j,l)
    -
    138  qm = amax1(qm,h1m12)
    -
    139  qmcvg= qcnvg(i,j)
    -
    140 !
    +
    132 ! SET TARGET PRESSURES.
    +
    133 
    +
    134  p100 = pint(i,j,nint(lmh(i,j)))
    +
    135  p98 = sig98*p100
    +
    136  p96 = sig96*p100
    +
    137  p85 = sig85*p100
    +
    138  p84 = sig84*p100
    +
    139  p47 = sig47*p100
    +
    140  p18 = sig18*p100
    141 !
    -
    142 ! COMPUTE RELATIVE HUMIDITY.
    -
    143 !
    -
    144  qs=pq0/pm*exp(a2*(tm-a3)/(tm-a4))
    -
    145 !
    -
    146  rh = qm/qs
    -
    147  IF (rh>h1) THEN
    -
    148  rh = h1
    -
    149  qm = rh*qs
    -
    150  ENDIF
    -
    151  IF (rh<d01) THEN
    -
    152  rh = d01
    -
    153  qm = rh*qs
    -
    154  ENDIF
    -
    155 !
    -
    156 ! SIGMA 0.85-1.00 MOISTURE CONVERGENCE.
    -
    157  IF ((pm<=p100).AND.(pm>=p85)) THEN
    -
    158  z8510(i,j) = z8510(i,j) + dz
    -
    159  qm8510(i,j) = qm8510(i,j) + qmcvg*dz
    -
    160  ENDIF
    -
    161 !
    -
    162 ! SIGMA 0.47-1.00 RELATIVE HUMIDITY.
    -
    163  IF ((pm<=p100).AND.(pm>=p47)) THEN
    -
    164  z4710(i,j) = z4710(i,j) + dz
    -
    165  rh4710(i,j) = rh4710(i,j) + rh*dz
    -
    166  ENDIF
    -
    167 !
    -
    168 ! SIGMA 0.84-0.98 RELATIVE HUMIDITY.
    -
    169  IF ((pm<=p98).AND.(pm>=p84)) THEN
    -
    170  z8498(i,j) = z8498(i,j) + dz
    -
    171  rh8498(i,j) = rh8498(i,j) + rh*dz
    -
    172  ENDIF
    -
    173 !
    -
    174 ! SIGMA 0.47-0.96 RELATIVE HUMIDITY.
    -
    175  IF ((pm<=p96).AND.(pm>=p47)) THEN
    -
    176  z4796(i,j) = z4796(i,j) + dz
    -
    177  rh4796(i,j) = rh4796(i,j) + rh*dz
    -
    178  ENDIF
    -
    179 !
    -
    180 ! SIGMA 0.18-0.47 RELATIVE HUMIDITY.
    -
    181  IF ((pm<=p47).AND.(pm>=p18)) THEN
    -
    182  z1847(i,j) = z1847(i,j) + dz
    -
    183  rh1847(i,j) = rh1847(i,j) + rh*dz
    -
    184  ENDIF
    -
    185 !
    -
    186  ENDDO
    -
    187  ENDDO
    -
    188  ENDDO
    -
    189 !
    -
    190  DO j=jsta_m2,jend_m2
    -
    191  DO i=ista_m,iend_m
    -
    192 ! NORMALIZE TO GET LAYER MEAN VALUES.
    -
    193  IF (z8510(i,j)>0) THEN
    -
    194  qm8510(i,j) = qm8510(i,j)/z8510(i,j)
    -
    195  ELSE
    -
    196  qm8510(i,j) = spval
    -
    197  ENDIF
    -
    198  IF (abs(qm8510(i,j)-spval)<small)qm8510(i,j)=h1m12
    -
    199 !
    -
    200  IF (z4710(i,j)>0) THEN
    -
    201  rh4710(i,j) = rh4710(i,j)/z4710(i,j)
    -
    202  ELSE
    -
    203  rh4710(i,j) = spval
    -
    204  ENDIF
    -
    205 !
    -
    206  IF (z8498(i,j)>0) THEN
    -
    207  rh8498(i,j) = rh8498(i,j)/z8498(i,j)
    +
    142 !
    +
    143 ! COMPUTE LAYER MEAN FIELDS AT THE GIVEN K.
    +
    144 !
    +
    145 ! COMPUTE P, Z, T, AND Q AT THE MIDPOINT OF THE CURRENT ETA LAYER.
    +
    146  alpm = d50*(alpint(i,j,l)+alpint(i,j,l+1))
    +
    147  dz = zint(i,j,l)-zint(i,j,l+1)
    +
    148  pm = exp(alpm)
    +
    149  tm = t(i,j,l)
    +
    150  qm = q(i,j,l)
    +
    151  qm = amax1(qm,h1m12)
    +
    152  qmcvg= qcnvg(i,j)
    +
    153 !
    +
    154 !
    +
    155 ! COMPUTE RELATIVE HUMIDITY.
    +
    156 !
    +
    157  qs=pq0/pm*exp(a2*(tm-a3)/(tm-a4))
    +
    158 !
    +
    159  rh = qm/qs
    +
    160  IF (rh>h1) THEN
    +
    161  rh = h1
    +
    162  qm = rh*qs
    +
    163  ENDIF
    +
    164  IF (rh<d01) THEN
    +
    165  rh = d01
    +
    166  qm = rh*qs
    +
    167  ENDIF
    +
    168 !
    +
    169 ! SIGMA 0.85-1.00 MOISTURE CONVERGENCE.
    +
    170  IF ((pm<=p100).AND.(pm>=p85)) THEN
    +
    171  z8510(i,j) = z8510(i,j) + dz
    +
    172  qm8510(i,j) = qm8510(i,j) + qmcvg*dz
    +
    173  ENDIF
    +
    174 !
    +
    175 ! SIGMA 0.47-1.00 RELATIVE HUMIDITY.
    +
    176  IF ((pm<=p100).AND.(pm>=p47)) THEN
    +
    177  z4710(i,j) = z4710(i,j) + dz
    +
    178  rh4710(i,j) = rh4710(i,j) + rh*dz
    +
    179  ENDIF
    +
    180 !
    +
    181 ! SIGMA 0.84-0.98 RELATIVE HUMIDITY.
    +
    182  IF ((pm<=p98).AND.(pm>=p84)) THEN
    +
    183  z8498(i,j) = z8498(i,j) + dz
    +
    184  rh8498(i,j) = rh8498(i,j) + rh*dz
    +
    185  ENDIF
    +
    186 !
    +
    187 ! SIGMA 0.47-0.96 RELATIVE HUMIDITY.
    +
    188  IF ((pm<=p96).AND.(pm>=p47)) THEN
    +
    189  z4796(i,j) = z4796(i,j) + dz
    +
    190  rh4796(i,j) = rh4796(i,j) + rh*dz
    +
    191  ENDIF
    +
    192 !
    +
    193 ! SIGMA 0.18-0.47 RELATIVE HUMIDITY.
    +
    194  IF ((pm<=p47).AND.(pm>=p18)) THEN
    +
    195  z1847(i,j) = z1847(i,j) + dz
    +
    196  rh1847(i,j) = rh1847(i,j) + rh*dz
    +
    197  ENDIF
    +
    198 !
    +
    199  ENDDO
    +
    200  ENDDO
    +
    201  ENDDO
    +
    202 !
    +
    203  DO j=jsta_m2,jend_m2
    +
    204  DO i=ista_m,iend_m
    +
    205 ! NORMALIZE TO GET LAYER MEAN VALUES.
    +
    206  IF (z8510(i,j)>0) THEN
    +
    207  qm8510(i,j) = qm8510(i,j)/z8510(i,j)
    208  ELSE
    -
    209  rh8498(i,j) = spval
    +
    209  qm8510(i,j) = spval
    210  ENDIF
    -
    211 !
    -
    212  IF (z4796(i,j)>0) THEN
    -
    213  rh4796(i,j) = rh4796(i,j)/z4796(i,j)
    -
    214  ELSE
    -
    215  rh4796(i,j) = spval
    -
    216  ENDIF
    -
    217 !
    -
    218  IF (z1847(i,j)>0) THEN
    -
    219  rh1847(i,j) = rh1847(i,j)/z1847(i,j)
    -
    220  ELSE
    -
    221  rh1847(i,j) = spval
    -
    222  ENDIF
    -
    223  ENDDO
    -
    224  ENDDO
    -
    225 !
    -
    226 !
    -
    227 ! END OF ROUTINE.
    -
    228 !
    -
    229  RETURN
    -
    230  END
    -
    231 
    +
    211  IF (abs(qm8510(i,j)-spval)<small)qm8510(i,j)=h1m12
    +
    212 !
    +
    213  IF (z4710(i,j)>0) THEN
    +
    214  rh4710(i,j) = rh4710(i,j)/z4710(i,j)
    +
    215  ELSE
    +
    216  rh4710(i,j) = spval
    +
    217  ENDIF
    +
    218 !
    +
    219  IF (z8498(i,j)>0) THEN
    +
    220  rh8498(i,j) = rh8498(i,j)/z8498(i,j)
    +
    221  ELSE
    +
    222  rh8498(i,j) = spval
    +
    223  ENDIF
    +
    224 !
    +
    225  IF (z4796(i,j)>0) THEN
    +
    226  rh4796(i,j) = rh4796(i,j)/z4796(i,j)
    +
    227  ELSE
    +
    228  rh4796(i,j) = spval
    +
    229  ENDIF
    +
    230 !
    +
    231  IF (z1847(i,j)>0) THEN
    +
    232  rh1847(i,j) = rh1847(i,j)/z1847(i,j)
    +
    233  ELSE
    +
    234  rh1847(i,j) = spval
    +
    235  ENDIF
    +
    236  ENDDO
    +
    237  ENDDO
    +
    238 !
    +
    239 !
    +
    240 ! END OF ROUTINE.
    +
    241 !
    +
    242  RETURN
    +
    243  END
    +
    244 
    +
    ctlblk_mod
    Definition: CTLBLK.f:1
    masks
    Definition: MASKS_mod.f:1
    params_mod
    Definition: params.F:1
    - + +
    ngmfld
    subroutine ngmfld(RH4710, RH4796, RH1847, RH8498, QM8510)
    ngmfld() computes layer mean NGM fields
    Definition: NGMFLD.f:59
    @@ -304,7 +309,7 @@
  • sorc
  • ncep_post.fd
  • NGMFLD.f
  • Generated by - doxygen 1.8.6
    • + doxygen 1.8.5 diff --git a/OTLFT_8f.html b/OTLFT_8f.html index e317ebb58..287f48749 100644 --- a/OTLFT_8f.html +++ b/OTLFT_8f.html @@ -3,7 +3,7 @@ - + UPP: OTLFT.f File Reference @@ -30,7 +30,7 @@
    UPP -  V11.0.0 +  11.0.0
    @@ -38,7 +38,7 @@ - + @@ -109,20 +109,20 @@
    -

    otlft() computes lifted index. +

    otlft() computes lifted index. More...

    Go to the source code of this file.

    - - + + +

    Functions/Subroutines

    -subroutine OTLFT (PBND, TBND, QBND, SLINDX)
     
    subroutine otlft (PBND, TBND, QBND, SLINDX)
     otlft() computes lifted index. More...
     

    Detailed Description

    -

    otlft() computes lifted index.

    -

    This routine computes lifts a parcel specified by the passed pressure, temperature, and specific humidity to 500mb and then computes a lifted index. This lifted lifted index is the difference between the lifted parcel's temperature at 500mb and the ambient 500mb temperature.

    +

    otlft() computes lifted index.

    +

    This routine lifts a parcel specified by the passed pressure, temperature, and specific humidity to 500mb and then computes a lifted index. This lifted index is the difference between the lifted parcel's temperature at 500mb and the ambient 500mb temperature.

    Parameters
    @@ -155,7 +155,64 @@

    Program History Log

    Date
    1993-03-10

    Definition in file OTLFT.f.

    - +

    Function/Subroutine Documentation

    + +
    +
    +
    [in]PBNDParcel pressure.
    + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
    subroutine otlft (real, dimension(ista:iend,jsta:jend), intent(in) PBND,
    real, dimension(ista:iend,jsta:jend), intent(in) TBND,
    real, dimension(ista:iend,jsta:jend), intent(in) QBND,
    real, dimension(ista:iend,jsta:jend), intent(out) SLINDX 
    )
    +
    + +

    otlft() computes lifted index.

    +

    This routine lifts a parcel specified by the passed pressure, temperature, and specific humidity to 500mb and then computes a lifted index. This lifted index is the difference between the lifted parcel's temperature at 500mb and the ambient 500mb temperature.

    +
    Parameters
    + + + + + +
    [in]PBNDParcel pressure.
    [in]TBNDParcel temperature.
    [in]QBNDParcel specific humidity.
    [out]SLINDXLifted index.
    +
    +
    + +

    Definition at line 44 of file OTLFT.f.

    + +

    References upp_physics::fpvsnew().

    + +

    Referenced by miscln().

    + +
    +
    +
    diff --git a/OTLFT_8f.js b/OTLFT_8f.js index 13db495cb..5638a5a17 100644 --- a/OTLFT_8f.js +++ b/OTLFT_8f.js @@ -1,4 +1,4 @@ var OTLFT_8f = [ - [ "OTLFT", "OTLFT_8f.html#a544edbe108f8617b421567a818a2fcbb", null ] + [ "otlft", "OTLFT_8f.html#a44baec1992cd116c348393de3daf2f2b", null ] ]; \ No newline at end of file diff --git a/OTLFT_8f_source.html b/OTLFT_8f_source.html index 0555b292b..a8107dfbc 100644 --- a/OTLFT_8f_source.html +++ b/OTLFT_8f_source.html @@ -3,7 +3,7 @@ - + UPP: OTLFT.f Source File @@ -30,7 +30,7 @@
    UPP -  V11.0.0 +  11.0.0
    @@ -38,7 +38,7 @@ - + @@ -107,210 +107,215 @@
    Go to the documentation of this file.
    1 
    -
    28  SUBROUTINE otlft(PBND,TBND,QBND,SLINDX)
    -
    29 
    -
    30 !
    -
    31 !
    -
    32  use vrbls2d, only: t500
    -
    33  use lookup_mod, only: thl, rdth, jtb, qs0, sqs, rdq, itb, ptbl, &
    -
    34  pl, rdp, the0, sthe, rdthe, ttbl
    -
    35  use ctlblk_mod, only: jsta, jend, im, spval, ista, iend
    -
    36  use params_mod, only: d00, h10e5, capa, elocp, eps, oneps
    -
    37  use upp_physics, only: fpvsnew
    -
    38 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    -
    39  implicit none
    -
    40 !
    -
    41 ! SET LOCAL PARAMETERS.
    -
    42  real,PARAMETER :: d8202=.820231e0 , h5e4=5.e4 , p500=50000.
    +
    28 
    +
    29 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    43 
    -
    44 !
    -
    45 ! DECLARE VARIABLES.
    -
    46  real,dimension(ista:iend,jsta:jend),intent(in) :: pbnd,tbnd,qbnd
    -
    47  real,dimension(ista:iend,jsta:jend),intent(out) :: slindx
    -
    48  REAL :: tvp, esatp, qsatp
    -
    49  REAL :: bqs00, sqs00, bqs10, sqs10, p00, p10, p01, p11, bq, sq, tq
    -
    50  REAL :: bthe00, sthe00, bthe10, sthe10, bth, sth, tth
    -
    51  REAL :: t00, t10, t01, t11, tbt, qbt, apebt, tthbt, ppq, pp
    -
    52  REAL :: tqq, qq, tpsp, apesp, tthes, tp, partmp
    -
    53 !
    -
    54  INTEGER :: i, j, ittbk, iq, it, iptbk, ith, ip
    -
    55  INTEGER :: ittb, iqtb, iptb, ithtb
    -
    56 !
    -
    57 !********************************************************************
    -
    58 ! START OTLFT HERE.
    -
    59 !
    -
    60 ! ZERO LIFTED INDEX ARRAY.
    -
    61 !
    -
    62 !$omp parallel do private(i,j)
    -
    63  DO j=jsta,jend
    -
    64  DO i=ista,iend
    -
    65  slindx(i,j) = d00
    -
    66  ENDDO
    -
    67  ENDDO
    -
    68 !
    -
    69 !--------------FIND EXNER IN BOUNDARY LAYER-----------------------------
    -
    70 !
    -
    71  DO j=jsta,jend
    -
    72  DO i=ista,iend
    -
    73  tbt = tbnd(i,j)
    -
    74  qbt = qbnd(i,j)
    -
    75 !
    -
    76  if( tbt < spval ) then
    -
    77 
    -
    78  apebt = (h10e5/pbnd(i,j))**capa
    -
    79 !
    -
    80 !--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX--------------
    -
    81 !
    -
    82  tthbt = tbt*apebt
    -
    83  tth=(tthbt-thl)*rdth
    -
    84  tqq = tth-aint(tth)
    -
    85  ittb = int(tth)+1
    +
    44  SUBROUTINE otlft(PBND,TBND,QBND,SLINDX)
    +
    45 
    +
    46 !
    +
    47 !
    +
    48  use vrbls2d, only: t500
    +
    49  use lookup_mod, only: thl, rdth, jtb, qs0, sqs, rdq, itb, ptbl, &
    +
    50  pl, rdp, the0, sthe, rdthe, ttbl
    +
    51  use ctlblk_mod, only: jsta, jend, im, spval, ista, iend
    +
    52  use params_mod, only: d00, h10e5, capa, elocp, eps, oneps
    +
    53  use upp_physics, only: fpvsnew
    +
    54 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    +
    55  implicit none
    +
    56 !
    +
    57 ! SET LOCAL PARAMETERS.
    +
    58  real,PARAMETER :: d8202=.820231e0 , h5e4=5.e4 , p500=50000.
    +
    59 
    +
    60 !
    +
    61 ! DECLARE VARIABLES.
    +
    62  real,dimension(ista:iend,jsta:jend),intent(in) :: pbnd,tbnd,qbnd
    +
    63  real,dimension(ista:iend,jsta:jend),intent(out) :: slindx
    +
    64  REAL :: tvp, esatp, qsatp
    +
    65  REAL :: bqs00, sqs00, bqs10, sqs10, p00, p10, p01, p11, bq, sq, tq
    +
    66  REAL :: bthe00, sthe00, bthe10, sthe10, bth, sth, tth
    +
    67  REAL :: t00, t10, t01, t11, tbt, qbt, apebt, tthbt, ppq, pp
    +
    68  REAL :: tqq, qq, tpsp, apesp, tthes, tp, partmp
    +
    69 !
    +
    70  INTEGER :: i, j, ittbk, iq, it, iptbk, ith, ip
    +
    71  INTEGER :: ittb, iqtb, iptb, ithtb
    +
    72 !
    +
    73 !********************************************************************
    +
    74 ! START OTLFT HERE.
    +
    75 !
    +
    76 ! ZERO LIFTED INDEX ARRAY.
    +
    77 !
    +
    78 !$omp parallel do private(i,j)
    +
    79  DO j=jsta,jend
    +
    80  DO i=ista,iend
    +
    81  slindx(i,j) = d00
    +
    82  ENDDO
    +
    83  ENDDO
    +
    84 !
    +
    85 !--------------FIND EXNER IN BOUNDARY LAYER-----------------------------
    86 !
    -
    87 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    -
    88 !
    -
    89  IF(ittb < 1)THEN
    -
    90  ittb = 1
    -
    91  tqq = d00
    -
    92  ENDIF
    -
    93  IF(ittb >= jtb)THEN
    -
    94  ittb = jtb-1
    -
    95  tqq = d00
    -
    96  ENDIF
    +
    87  DO j=jsta,jend
    +
    88  DO i=ista,iend
    +
    89  tbt = tbnd(i,j)
    +
    90  qbt = qbnd(i,j)
    +
    91 !
    +
    92  if( tbt < spval ) then
    +
    93 
    +
    94  apebt = (h10e5/pbnd(i,j))**capa
    +
    95 !
    +
    96 !--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX--------------
    97 !
    -
    98 !--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
    -
    99 !
    -
    100  ittbk = ittb
    -
    101  bqs00=qs0(ittbk)
    -
    102  sqs00=sqs(ittbk)
    -
    103  bqs10=qs0(ittbk+1)
    -
    104  sqs10=sqs(ittbk+1)
    -
    105 !
    -
    106 !--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
    -
    107 !
    -
    108  bq=(bqs10-bqs00)*tqq+bqs00
    -
    109  sq=(sqs10-sqs00)*tqq+sqs00
    -
    110  tq=(qbt-bq)/sq*rdq
    -
    111  ppq = tq-aint(tq)
    -
    112  iqtb = int(tq)+1
    +
    98  tthbt = tbt*apebt
    +
    99  tth=(tthbt-thl)*rdth
    +
    100  tqq = tth-aint(tth)
    +
    101  ittb = int(tth)+1
    +
    102 !
    +
    103 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    +
    104 !
    +
    105  IF(ittb < 1)THEN
    +
    106  ittb = 1
    +
    107  tqq = d00
    +
    108  ENDIF
    +
    109  IF(ittb >= jtb)THEN
    +
    110  ittb = jtb-1
    +
    111  tqq = d00
    +
    112  ENDIF
    113 !
    -
    114 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    +
    114 !--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
    115 !
    -
    116  IF(iqtb < 1)THEN
    -
    117  iqtb = 1
    -
    118  ppq = d00
    -
    119  ENDIF
    -
    120  IF(iqtb >= itb)THEN
    -
    121  iqtb = itb-1
    -
    122  ppq = d00
    -
    123  ENDIF
    -
    124 !
    -
    125 !--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
    -
    126 !
    -
    127  iq=iqtb
    -
    128  it=ittb
    -
    129  p00=ptbl(iq,it)
    -
    130  p10=ptbl(iq+1,it)
    -
    131  p01=ptbl(iq,it+1)
    -
    132  p11=ptbl(iq+1,it+1)
    -
    133 !
    -
    134 !--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
    -
    135 !
    -
    136  tpsp = p00+(p10-p00)*ppq+(p01-p00)*tqq &
    -
    137  +(p00-p10-p01+p11)*ppq*tqq
    -
    138  IF(tpsp <= d00) tpsp = h10e5
    -
    139  apesp = (h10e5/tpsp)**capa
    -
    140  tthes = tthbt*exp(elocp*qbt*apesp/tthbt)
    -
    141 !
    -
    142 !-----------------------------------------------------------------------
    -
    143 !
    -
    144 !
    -
    145 !--------------SCALING PRESSURE & TT TABLE INDEX------------------------
    -
    146 !
    -
    147  tp = (h5e4-pl)*rdp
    -
    148  qq = tp-aint(tp)
    -
    149  iptb = int(tp)+1
    -
    150 !
    -
    151 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    -
    152 !
    -
    153  IF(iptb < 1)THEN
    -
    154  iptb = 1
    -
    155  qq = d00
    -
    156  ENDIF
    -
    157  IF(iptb >= itb)THEN
    -
    158  iptb = itb-1
    -
    159  qq = d00
    -
    160  ENDIF
    -
    161 !
    -
    162 !--------------BASE AND SCALING FACTOR FOR THE--------------------------
    -
    163 !
    -
    164  iptbk=iptb
    -
    165  bthe00=the0(iptbk)
    -
    166  sthe00=sthe(iptbk)
    -
    167  bthe10=the0(iptbk+1)
    -
    168  sthe10=sthe(iptbk+1)
    -
    169 !
    -
    170 !--------------SCALING THE & TT TABLE INDEX-----------------------------
    -
    171 !
    -
    172  bth=(bthe10-bthe00)*qq+bthe00
    -
    173  sth=(sthe10-sthe00)*qq+sthe00
    -
    174  tth=(tthes-bth)/sth*rdthe
    -
    175  pp = tth-aint(tth)
    -
    176  ithtb = int(tth)+1
    +
    116  ittbk = ittb
    +
    117  bqs00=qs0(ittbk)
    +
    118  sqs00=sqs(ittbk)
    +
    119  bqs10=qs0(ittbk+1)
    +
    120  sqs10=sqs(ittbk+1)
    +
    121 !
    +
    122 !--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
    +
    123 !
    +
    124  bq=(bqs10-bqs00)*tqq+bqs00
    +
    125  sq=(sqs10-sqs00)*tqq+sqs00
    +
    126  tq=(qbt-bq)/sq*rdq
    +
    127  ppq = tq-aint(tq)
    +
    128  iqtb = int(tq)+1
    +
    129 !
    +
    130 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    +
    131 !
    +
    132  IF(iqtb < 1)THEN
    +
    133  iqtb = 1
    +
    134  ppq = d00
    +
    135  ENDIF
    +
    136  IF(iqtb >= itb)THEN
    +
    137  iqtb = itb-1
    +
    138  ppq = d00
    +
    139  ENDIF
    +
    140 !
    +
    141 !--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
    +
    142 !
    +
    143  iq=iqtb
    +
    144  it=ittb
    +
    145  p00=ptbl(iq,it)
    +
    146  p10=ptbl(iq+1,it)
    +
    147  p01=ptbl(iq,it+1)
    +
    148  p11=ptbl(iq+1,it+1)
    +
    149 !
    +
    150 !--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
    +
    151 !
    +
    152  tpsp = p00+(p10-p00)*ppq+(p01-p00)*tqq &
    +
    153  +(p00-p10-p01+p11)*ppq*tqq
    +
    154  IF(tpsp <= d00) tpsp = h10e5
    +
    155  apesp = (h10e5/tpsp)**capa
    +
    156  tthes = tthbt*exp(elocp*qbt*apesp/tthbt)
    +
    157 !
    +
    158 !-----------------------------------------------------------------------
    +
    159 !
    +
    160 !
    +
    161 !--------------SCALING PRESSURE & TT TABLE INDEX------------------------
    +
    162 !
    +
    163  tp = (h5e4-pl)*rdp
    +
    164  qq = tp-aint(tp)
    +
    165  iptb = int(tp)+1
    +
    166 !
    +
    167 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    +
    168 !
    +
    169  IF(iptb < 1)THEN
    +
    170  iptb = 1
    +
    171  qq = d00
    +
    172  ENDIF
    +
    173  IF(iptb >= itb)THEN
    +
    174  iptb = itb-1
    +
    175  qq = d00
    +
    176  ENDIF
    177 !
    -
    178 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    +
    178 !--------------BASE AND SCALING FACTOR FOR THE--------------------------
    179 !
    -
    180  IF(ithtb < 1)THEN
    -
    181  ithtb = 1
    -
    182  pp = d00
    -
    183  ENDIF
    -
    184  IF(ithtb >= jtb)THEN
    -
    185  ithtb = jtb-1
    -
    186  pp = d00
    -
    187  ENDIF
    -
    188 !
    -
    189 !--------------TEMPERATURE AT FOUR SURROUNDING TT TABLE PTS.------------
    -
    190 !
    -
    191  ith=ithtb
    -
    192  ip=iptb
    -
    193  t00=ttbl(ith,ip)
    -
    194  t10=ttbl(ith+1,ip)
    -
    195  t01=ttbl(ith,ip+1)
    -
    196  t11=ttbl(ith+1,ip+1)
    -
    197 !
    -
    198 !--------------PARCEL TEMPERATURE AT 500MB----------------------------
    -
    199 !
    -
    200  IF(tpsp >= h5e4)THEN
    -
    201  partmp=(t00+(t10-t00)*pp+(t01-t00)*qq &
    -
    202  +(t00-t10-t01+t11)*pp*qq)
    -
    203  ELSE
    -
    204  partmp=tbt*apebt*d8202
    -
    205  ENDIF
    +
    180  iptbk=iptb
    +
    181  bthe00=the0(iptbk)
    +
    182  sthe00=sthe(iptbk)
    +
    183  bthe10=the0(iptbk+1)
    +
    184  sthe10=sthe(iptbk+1)
    +
    185 !
    +
    186 !--------------SCALING THE & TT TABLE INDEX-----------------------------
    +
    187 !
    +
    188  bth=(bthe10-bthe00)*qq+bthe00
    +
    189  sth=(sthe10-sthe00)*qq+sthe00
    +
    190  tth=(tthes-bth)/sth*rdthe
    +
    191  pp = tth-aint(tth)
    +
    192  ithtb = int(tth)+1
    +
    193 !
    +
    194 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    +
    195 !
    +
    196  IF(ithtb < 1)THEN
    +
    197  ithtb = 1
    +
    198  pp = d00
    +
    199  ENDIF
    +
    200  IF(ithtb >= jtb)THEN
    +
    201  ithtb = jtb-1
    +
    202  pp = d00
    +
    203  ENDIF
    +
    204 !
    +
    205 !--------------TEMPERATURE AT FOUR SURROUNDING TT TABLE PTS.------------
    206 !
    -
    207 !--------------LIFTED INDEX---------------------------------------------
    -
    208 !
    -
    209 ! GSM THE PARCEL TEMPERATURE AT 500 MB HAS BEEN COMPUTED, AND WE
    -
    210 ! FIND THE MIXING RATIO AT THAT LEVEL WHICH WILL BE THE SATURATION
    -
    211 ! VALUE SINCE WE'RE FOLLOWING A MOIST ADIABAT. NOTE THAT THE
    -
    212 ! AMBIENT 500 MB SHOULD PROBABLY BE VIRTUALIZED, BUT THE IMPACT
    -
    213 ! OF MOISTURE AT THAT LEVEL IS QUITE SMALL
    -
    214  esatp=fpvsnew(partmp)
    -
    215  qsatp=eps*esatp/(p500-esatp*oneps)
    -
    216  tvp=partmp*(1+0.608*qsatp)
    -
    217  slindx(i,j)=t500(i,j)-tvp
    -
    218 
    -
    219  else
    -
    220  slindx(i,j)=spval
    -
    221  endif
    -
    222  END DO
    -
    223  END DO
    -
    224 !
    -
    225 ! END OF ROUTINE.
    -
    226  RETURN
    -
    227  END
    +
    207  ith=ithtb
    +
    208  ip=iptb
    +
    209  t00=ttbl(ith,ip)
    +
    210  t10=ttbl(ith+1,ip)
    +
    211  t01=ttbl(ith,ip+1)
    +
    212  t11=ttbl(ith+1,ip+1)
    +
    213 !
    +
    214 !--------------PARCEL TEMPERATURE AT 500MB----------------------------
    +
    215 !
    +
    216  IF(tpsp >= h5e4)THEN
    +
    217  partmp=(t00+(t10-t00)*pp+(t01-t00)*qq &
    +
    218  +(t00-t10-t01+t11)*pp*qq)
    +
    219  ELSE
    +
    220  partmp=tbt*apebt*d8202
    +
    221  ENDIF
    +
    222 !
    +
    223 !--------------LIFTED INDEX---------------------------------------------
    +
    224 !
    +
    225 ! GSM THE PARCEL TEMPERATURE AT 500 MB HAS BEEN COMPUTED, AND WE
    +
    226 ! FIND THE MIXING RATIO AT THAT LEVEL WHICH WILL BE THE SATURATION
    +
    227 ! VALUE SINCE WE'RE FOLLOWING A MOIST ADIABAT. NOTE THAT THE
    +
    228 ! AMBIENT 500 MB SHOULD PROBABLY BE VIRTUALIZED, BUT THE IMPACT
    +
    229 ! OF MOISTURE AT THAT LEVEL IS QUITE SMALL
    +
    230  esatp=fpvsnew(partmp)
    +
    231  qsatp=eps*esatp/(p500-esatp*oneps)
    +
    232  tvp=partmp*(1+0.608*qsatp)
    +
    233  slindx(i,j)=t500(i,j)-tvp
    +
    234 
    +
    235  else
    +
    236  slindx(i,j)=spval
    +
    237  endif
    +
    238  END DO
    +
    239  END DO
    +
    240 !
    +
    241 ! END OF ROUTINE.
    +
    242  RETURN
    +
    243  END
    +
    ctlblk_mod
    Definition: CTLBLK.f:1
    params_mod
    Definition: params.F:1
    -
    upp_physics::fpvsnew
    elemental real function, public fpvsnew(t)
    Definition: UPP_PHYSICS.f:345
    -
    upp_physics
    calcape() computes CAPE/CINS and other storm related variables.
    Definition: UPP_PHYSICS.f:27
    +
    otlft
    subroutine otlft(PBND, TBND, QBND, SLINDX)
    otlft() computes lifted index.
    Definition: OTLFT.f:44
    +
    upp_physics::fpvsnew
    elemental real function, public fpvsnew(t)
    Definition: UPP_PHYSICS.f:359
    +
    upp_physics
    Definition: UPP_PHYSICS.f:40
    lookup_mod
    Definition: LOOKUP.f:1
    @@ -320,7 +325,7 @@
  • sorc
  • ncep_post.fd
  • OTLFT.f
  • Generated by - doxygen 1.8.6
    • + doxygen 1.8.5 diff --git a/OTLIFT_8f.html b/OTLIFT_8f.html index 7309d0f6c..b0c4f0555 100644 --- a/OTLIFT_8f.html +++ b/OTLIFT_8f.html @@ -3,7 +3,7 @@ - + UPP: OTLIFT.f File Reference @@ -30,7 +30,7 @@
    UPP -  V11.0.0 +  11.0.0
    @@ -38,7 +38,7 @@ - + @@ -109,20 +109,20 @@
    -

    otlift() computes SFC to 500mb lifted index. +

    otlift() computes SFC to 500mb lifted index. More...

    Go to the source code of this file.

    - - + + +

    Functions/Subroutines

    -subroutine OTLIFT (SLINDX)
     
    subroutine otlift (SLINDX)
     otlift() computes SFC to 500mb lifted index. More...
     

    Detailed Description

    -

    otlift() computes SFC to 500mb lifted index.

    -

    This routine computes a surface to 500mb lifted index. The lifted parcel is from the first atmpspheric ETA layer (ie, the ETA layer closest to the model ground). The lifted index is the difference between this parcel's temperature at 500mb and the ambient 500mb temperature.

    +

    otlift() computes SFC to 500mb lifted index.

    +

    This routine computes a surface to 500mb lifted index. The lifted parcel is from the first atmospheric ETA layer (ie, the ETA layer closest to the model ground). The lifted index is the difference between this parcel's temperature at 500mb and the ambient 500mb temperature.

    Parameters
    @@ -156,7 +156,37 @@

    Program History Log

    Date
    1993-03-10

    Definition in file OTLIFT.f.

    - +

    Function/Subroutine Documentation

    + +
    +
    +
    [out]SLINDXlifted index.
    + + + + + + + +
    subroutine otlift (real, dimension(ista:iend,jsta:jend), intent(out) SLINDX)
    +
    + +

    otlift() computes SFC to 500mb lifted index.

    +

    This routine computes a surface to 500mb lifted index. The lifted parcel is from the first atmospheric ETA layer (ie, the ETA layer closest to the model ground). The lifted index is the difference between this parcel's temperature at 500mb and the ambient 500mb temperature.

    +
    Parameters
    + + +
    [out]SLINDXlifted index.
    +
    +
    + +

    Definition at line 38 of file OTLIFT.f.

    + +

    References upp_physics::fpvsnew().

    + +
    +
    +
    diff --git a/OTLIFT_8f.js b/OTLIFT_8f.js index bbbd8632b..34d4648cd 100644 --- a/OTLIFT_8f.js +++ b/OTLIFT_8f.js @@ -1,4 +1,4 @@ var OTLIFT_8f = [ - [ "OTLIFT", "OTLIFT_8f.html#ac593a8e8f21ab04345f692158107e5f5", null ] + [ "otlift", "OTLIFT_8f.html#a5ad053b19cdb3e7a7ebc7a4a0d9642f3", null ] ]; \ No newline at end of file diff --git a/OTLIFT_8f_source.html b/OTLIFT_8f_source.html index 383b611e8..de3aa5aab 100644 --- a/OTLIFT_8f_source.html +++ b/OTLIFT_8f_source.html @@ -3,7 +3,7 @@ - + UPP: OTLIFT.f Source File @@ -30,7 +30,7 @@
    UPP -  V11.0.0 +  11.0.0
    @@ -38,7 +38,7 @@ - + @@ -107,178 +107,183 @@
    Go to the documentation of this file.
    1 
    -
    26  SUBROUTINE otlift(SLINDX)
    -
    27 
    -
    28 !
    -
    29  use vrbls3d, only: pmid, t, q
    -
    30  use vrbls2d, only: t500
    -
    31  use masks, only: lmh
    -
    32  use lookup_mod, only: thl, rdth, jtb, qs0, sqs, rdq,itb, ptbl, pl, &
    -
    33  rdp, the0, sthe, rdthe, ttbl
    -
    34  use ctlblk_mod, only: jsta, jend, im, spval, ista, iend
    -
    35  use params_mod, only: d00,h10e5, capa, elocp, eps, oneps
    -
    36  use upp_physics, only: fpvsnew
    -
    37 !
    -
    38 
    -
    39 !
    -
    40  implicit none
    -
    41 !
    -
    42 ! SET LOCAL PARAMETERS.
    -
    43  real,PARAMETER :: d8202=.820231e0 , h5e4=5.e4 , p500=50000.
    -
    44 
    -
    45 !
    -
    46 ! DECLARE VARIABLES.
    -
    47  real,intent(out) :: slindx(ista:iend,jsta:jend)
    -
    48  REAL :: tvp, esatp, qsatp
    -
    49  REAL :: tth, tp, apesp, partmp, thesp, tpsp
    -
    50  REAL :: bqs00, sqs00, bqs10, sqs10, bq, sq, tq
    -
    51  REAL :: p00, p10, p01, p11, t00, t10, t01, t11
    -
    52  REAL :: bthe00, sthe00, bthe10, sthe10, bth, sth
    -
    53  REAL :: tqq, qq, qbt, tthbt, tbt, apebt, ppq, pp
    -
    54 !
    -
    55  INTEGER :: i, j, lbtm, ittbk, iq, it, iptbk, ith, ip, iqtb
    -
    56  INTEGER :: ittb, iptb, ithtb
    +
    26 
    +
    27 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    +
    37 
    +
    38  SUBROUTINE otlift(SLINDX)
    +
    39 
    +
    40 !
    +
    41  use vrbls3d, only: pmid, t, q
    +
    42  use vrbls2d, only: t500
    +
    43  use masks, only: lmh
    +
    44  use lookup_mod, only: thl, rdth, jtb, qs0, sqs, rdq,itb, ptbl, pl, &
    +
    45  rdp, the0, sthe, rdthe, ttbl
    +
    46  use ctlblk_mod, only: jsta, jend, im, spval, ista, iend
    +
    47  use params_mod, only: d00,h10e5, capa, elocp, eps, oneps
    +
    48  use upp_physics, only: fpvsnew
    +
    49 !
    +
    50 
    +
    51 !
    +
    52  implicit none
    +
    53 !
    +
    54 ! SET LOCAL PARAMETERS.
    +
    55  real,PARAMETER :: d8202=.820231e0 , h5e4=5.e4 , p500=50000.
    +
    56 
    57 !
    -
    58 !***********************************************************************
    -
    59 ! START OTLIFT HERE
    -
    60 !
    -
    61 ! INTIALIZE LIFTED INDEX ARRAY TO ZERO.
    -
    62 !$omp parallel do private(i,j)
    -
    63  DO j=jsta,jend
    -
    64  DO i=ista,iend
    -
    65  slindx(i,j) = d00
    -
    66  ENDDO
    -
    67  ENDDO
    -
    68 !--------------FIND EXNER AT LOWEST LEVEL-------------------------------
    -
    69  DO j=jsta,jend
    -
    70  DO i=ista,iend
    -
    71  lbtm=nint(lmh(i,j))
    -
    72  IF(t(i,j,lbtm)<spval .AND. q(i,j,lbtm)<spval) THEN
    -
    73  tbt = t(i,j,lbtm)
    -
    74  qbt = q(i,j,lbtm)
    -
    75  apebt = (h10e5/pmid(i,j,lbtm))**capa
    -
    76 !--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX--------------
    -
    77  tthbt = tbt*apebt
    -
    78  tth = (tthbt-thl)*rdth
    -
    79  tqq = tth-aint(tth)
    -
    80  ittb = int(tth)+1
    -
    81 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    -
    82  IF(ittb < 1)THEN
    -
    83  ittb = 1
    -
    84  tqq = d00
    -
    85  ENDIF
    -
    86  IF(ittb >= jtb)THEN
    -
    87  ittb = jtb-1
    -
    88  tqq = d00
    -
    89  ENDIF
    -
    90 !--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
    -
    91  ittbk = ittb
    -
    92  bqs00=qs0(ittbk)
    -
    93  sqs00=sqs(ittbk)
    -
    94  bqs10=qs0(ittbk+1)
    -
    95  sqs10=sqs(ittbk+1)
    -
    96 !--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
    -
    97  bq=(bqs10-bqs00)*tqq+bqs00
    -
    98  sq=(sqs10-sqs00)*tqq+sqs00
    -
    99  tq=(qbt-bq)/sq*rdq
    -
    100  ppq = tq-aint(tq)
    -
    101  iqtb = int(tq)+1
    -
    102 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    -
    103  IF(iqtb < 1)THEN
    -
    104  iqtb = 1
    -
    105  ppq = d00
    -
    106  ENDIF
    -
    107  IF(iqtb >= itb)THEN
    -
    108  iqtb = itb-1
    -
    109  ppq = d00
    -
    110  ENDIF
    -
    111 !--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
    -
    112  iq=iqtb
    -
    113  it = ittb
    -
    114  p00=ptbl(iq,it)
    -
    115  p10=ptbl(iq+1,it)
    -
    116  p01=ptbl(iq,it+1)
    -
    117  p11=ptbl(iq+1,it+1)
    -
    118 !--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
    -
    119  tpsp = p00+(p10-p00)*ppq+(p01-p00)*tqq &
    -
    120  +(p00-p10-p01+p11)*ppq*tqq
    -
    121  IF(tpsp <= d00) tpsp = h10e5
    -
    122  apesp = (h10e5/tpsp)**capa
    -
    123  thesp = tthbt*exp(elocp*qbt*apesp/tthbt)
    -
    124 !--------------SCALING PRESSURE & TT TABLE INDEX------------------------
    -
    125  tp=(h5e4-pl)*rdp
    -
    126  qq = tp-aint(tp)
    -
    127  iptb = int(tp)+1
    -
    128 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    -
    129  IF(iptb < 1)THEN
    -
    130  iptb = 1
    -
    131  qq = d00
    -
    132  ENDIF
    -
    133  IF(iptb >= itb)THEN
    -
    134  iptb = itb-1
    -
    135  qq = d00
    -
    136  ENDIF
    -
    137 !--------------BASE AND SCALING FACTOR FOR THE--------------------------
    -
    138  iptbk=iptb
    -
    139  bthe00=the0(iptbk)
    -
    140  sthe00=sthe(iptbk)
    -
    141  bthe10=the0(iptbk+1)
    -
    142  sthe10=sthe(iptbk+1)
    -
    143 !--------------SCALING THE & TT TABLE INDEX-----------------------------
    -
    144  bth=(bthe10-bthe00)*qq+bthe00
    -
    145  sth=(sthe10-sthe00)*qq+sthe00
    -
    146  tth=(thesp-bth)/sth*rdthe
    -
    147  pp = tth-aint(tth)
    -
    148  ithtb = int(tth)+1
    -
    149 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    -
    150  IF(ithtb < 1)THEN
    -
    151  ithtb = 1
    -
    152  pp = d00
    -
    153  ENDIF
    -
    154  IF(ithtb >= jtb)THEN
    -
    155  ithtb = jtb-1
    -
    156  pp = d00
    -
    157  ENDIF
    -
    158 !--------------TEMPERATURE AT FOUR SURROUNDING TT TABLE PTS.------------
    -
    159  ith=ithtb
    -
    160  ip=iptb
    -
    161  t00=ttbl(ith,ip)
    -
    162  t10=ttbl(ith+1,ip)
    -
    163  t01=ttbl(ith,ip+1)
    -
    164  t11=ttbl(ith+1,ip+1)
    -
    165 !--------------PARCEL TEMPERATURE AT 500MB----------------------------
    -
    166  IF(tpsp >= h5e4)THEN
    -
    167  partmp=(t00+(t10-t00)*pp+(t01-t00)*qq &
    -
    168  +(t00-t10-t01+t11)*pp*qq)
    -
    169  ELSE
    -
    170  partmp=tbt*apebt*d8202
    -
    171  ENDIF
    -
    172 !--------------LIFTED INDEX---------------------------------------------
    -
    173 !
    -
    174 ! GSM THE PARCEL TEMPERATURE AT 500 MB HAS BEEN COMPUTED, AND WE
    -
    175 ! FIND THE MIXING RATIO AT THAT LEVEL WHICH WILL BE THE SATURATION
    -
    176 ! VALUE SINCE WE'RE FOLLOWING A MOIST ADIABAT. NOTE THAT THE
    -
    177 ! AMBIENT 500 MB SHOULD PROBABLY BE VIRTUALIZED, BUT THE IMPACT
    -
    178 ! OF MOISTURE AT THAT LEVEL IS QUITE SMALL
    -
    179 
    -
    180  esatp=fpvsnew(partmp)
    -
    181  qsatp=eps*esatp/(p500-esatp*oneps)
    -
    182  tvp=partmp*(1+0.608*qsatp)
    -
    183  slindx(i,j)=t500(i,j)-tvp
    -
    184  ENDIF !end T(I,J,LBTM)<SPVAL
    -
    185  ENDDO
    -
    186  ENDDO
    -
    187 ! write(*,*) ' in otlift t500 partmp ',t500(1,1),partmp(1,1)
    -
    188 ! write(*,*) ' in otlift tbt ',tbt(1,1)
    -
    189 !
    -
    190  RETURN
    -
    191  END
    +
    58 ! DECLARE VARIABLES.
    +
    59  real,intent(out) :: slindx(ista:iend,jsta:jend)
    +
    60  REAL :: tvp, esatp, qsatp
    +
    61  REAL :: tth, tp, apesp, partmp, thesp, tpsp
    +
    62  REAL :: bqs00, sqs00, bqs10, sqs10, bq, sq, tq
    +
    63  REAL :: p00, p10, p01, p11, t00, t10, t01, t11
    +
    64  REAL :: bthe00, sthe00, bthe10, sthe10, bth, sth
    +
    65  REAL :: tqq, qq, qbt, tthbt, tbt, apebt, ppq, pp
    +
    66 !
    +
    67  INTEGER :: i, j, lbtm, ittbk, iq, it, iptbk, ith, ip, iqtb
    +
    68  INTEGER :: ittb, iptb, ithtb
    +
    69 !
    +
    70 !***********************************************************************
    +
    71 ! START OTLIFT HERE
    +
    72 !
    +
    73 ! INTIALIZE LIFTED INDEX ARRAY TO ZERO.
    +
    74 !$omp parallel do private(i,j)
    +
    75  DO j=jsta,jend
    +
    76  DO i=ista,iend
    +
    77  slindx(i,j) = d00
    +
    78  ENDDO
    +
    79  ENDDO
    +
    80 !--------------FIND EXNER AT LOWEST LEVEL-------------------------------
    +
    81  DO j=jsta,jend
    +
    82  DO i=ista,iend
    +
    83  lbtm=nint(lmh(i,j))
    +
    84  IF(t(i,j,lbtm)<spval .AND. q(i,j,lbtm)<spval) THEN
    +
    85  tbt = t(i,j,lbtm)
    +
    86  qbt = q(i,j,lbtm)
    +
    87  apebt = (h10e5/pmid(i,j,lbtm))**capa
    +
    88 !--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX--------------
    +
    89  tthbt = tbt*apebt
    +
    90  tth = (tthbt-thl)*rdth
    +
    91  tqq = tth-aint(tth)
    +
    92  ittb = int(tth)+1
    +
    93 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    +
    94  IF(ittb < 1)THEN
    +
    95  ittb = 1
    +
    96  tqq = d00
    +
    97  ENDIF
    +
    98  IF(ittb >= jtb)THEN
    +
    99  ittb = jtb-1
    +
    100  tqq = d00
    +
    101  ENDIF
    +
    102 !--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
    +
    103  ittbk = ittb
    +
    104  bqs00=qs0(ittbk)
    +
    105  sqs00=sqs(ittbk)
    +
    106  bqs10=qs0(ittbk+1)
    +
    107  sqs10=sqs(ittbk+1)
    +
    108 !--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
    +
    109  bq=(bqs10-bqs00)*tqq+bqs00
    +
    110  sq=(sqs10-sqs00)*tqq+sqs00
    +
    111  tq=(qbt-bq)/sq*rdq
    +
    112  ppq = tq-aint(tq)
    +
    113  iqtb = int(tq)+1
    +
    114 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    +
    115  IF(iqtb < 1)THEN
    +
    116  iqtb = 1
    +
    117  ppq = d00
    +
    118  ENDIF
    +
    119  IF(iqtb >= itb)THEN
    +
    120  iqtb = itb-1
    +
    121  ppq = d00
    +
    122  ENDIF
    +
    123 !--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
    +
    124  iq=iqtb
    +
    125  it = ittb
    +
    126  p00=ptbl(iq,it)
    +
    127  p10=ptbl(iq+1,it)
    +
    128  p01=ptbl(iq,it+1)
    +
    129  p11=ptbl(iq+1,it+1)
    +
    130 !--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
    +
    131  tpsp = p00+(p10-p00)*ppq+(p01-p00)*tqq &
    +
    132  +(p00-p10-p01+p11)*ppq*tqq
    +
    133  IF(tpsp <= d00) tpsp = h10e5
    +
    134  apesp = (h10e5/tpsp)**capa
    +
    135  thesp = tthbt*exp(elocp*qbt*apesp/tthbt)
    +
    136 !--------------SCALING PRESSURE & TT TABLE INDEX------------------------
    +
    137  tp=(h5e4-pl)*rdp
    +
    138  qq = tp-aint(tp)
    +
    139  iptb = int(tp)+1
    +
    140 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    +
    141  IF(iptb < 1)THEN
    +
    142  iptb = 1
    +
    143  qq = d00
    +
    144  ENDIF
    +
    145  IF(iptb >= itb)THEN
    +
    146  iptb = itb-1
    +
    147  qq = d00
    +
    148  ENDIF
    +
    149 !--------------BASE AND SCALING FACTOR FOR THE--------------------------
    +
    150  iptbk=iptb
    +
    151  bthe00=the0(iptbk)
    +
    152  sthe00=sthe(iptbk)
    +
    153  bthe10=the0(iptbk+1)
    +
    154  sthe10=sthe(iptbk+1)
    +
    155 !--------------SCALING THE & TT TABLE INDEX-----------------------------
    +
    156  bth=(bthe10-bthe00)*qq+bthe00
    +
    157  sth=(sthe10-sthe00)*qq+sthe00
    +
    158  tth=(thesp-bth)/sth*rdthe
    +
    159  pp = tth-aint(tth)
    +
    160  ithtb = int(tth)+1
    +
    161 !--------------KEEPING INDICES WITHIN THE TABLE-------------------------
    +
    162  IF(ithtb < 1)THEN
    +
    163  ithtb = 1
    +
    164  pp = d00
    +
    165  ENDIF
    +
    166  IF(ithtb >= jtb)THEN
    +
    167  ithtb = jtb-1
    +
    168  pp = d00
    +
    169  ENDIF
    +
    170 !--------------TEMPERATURE AT FOUR SURROUNDING TT TABLE PTS.------------
    +
    171  ith=ithtb
    +
    172  ip=iptb
    +
    173  t00=ttbl(ith,ip)
    +
    174  t10=ttbl(ith+1,ip)
    +
    175  t01=ttbl(ith,ip+1)
    +
    176  t11=ttbl(ith+1,ip+1)
    +
    177 !--------------PARCEL TEMPERATURE AT 500MB----------------------------
    +
    178  IF(tpsp >= h5e4)THEN
    +
    179  partmp=(t00+(t10-t00)*pp+(t01-t00)*qq &
    +
    180  +(t00-t10-t01+t11)*pp*qq)
    +
    181  ELSE
    +
    182  partmp=tbt*apebt*d8202
    +
    183  ENDIF
    +
    184 !--------------LIFTED INDEX---------------------------------------------
    +
    185 !
    +
    186 ! GSM THE PARCEL TEMPERATURE AT 500 MB HAS BEEN COMPUTED, AND WE
    +
    187 ! FIND THE MIXING RATIO AT THAT LEVEL WHICH WILL BE THE SATURATION
    +
    188 ! VALUE SINCE WE'RE FOLLOWING A MOIST ADIABAT. NOTE THAT THE
    +
    189 ! AMBIENT 500 MB SHOULD PROBABLY BE VIRTUALIZED, BUT THE IMPACT
    +
    190 ! OF MOISTURE AT THAT LEVEL IS QUITE SMALL
    +
    191 
    +
    192  esatp=fpvsnew(partmp)
    +
    193  qsatp=eps*esatp/(p500-esatp*oneps)
    +
    194  tvp=partmp*(1+0.608*qsatp)
    +
    195  slindx(i,j)=t500(i,j)-tvp
    +
    196  ENDIF !end T(I,J,LBTM)<SPVAL
    +
    197  ENDDO
    +
    198  ENDDO
    +
    199 ! write(*,*) ' in otlift t500 partmp ',t500(1,1),partmp(1,1)
    +
    200 ! write(*,*) ' in otlift tbt ',tbt(1,1)
    +
    201 !
    +
    202  RETURN
    +
    203  END
    +
    ctlblk_mod
    Definition: CTLBLK.f:1
    masks
    Definition: MASKS_mod.f:1
    params_mod
    Definition: params.F:1
    - -
    upp_physics::fpvsnew
    elemental real function, public fpvsnew(t)
    Definition: UPP_PHYSICS.f:345
    -
    upp_physics
    calcape() computes CAPE/CINS and other storm related variables.
    Definition: UPP_PHYSICS.f:27
    +
    otlift
    subroutine otlift(SLINDX)
    otlift() computes SFC to 500mb lifted index.
    Definition: OTLIFT.f:38
    + +
    upp_physics::fpvsnew
    elemental real function, public fpvsnew(t)
    Definition: UPP_PHYSICS.f:359
    +
    upp_physics
    Definition: UPP_PHYSICS.f:40
    lookup_mod
    Definition: LOOKUP.f:1
    @@ -288,7 +293,7 @@
  • sorc
  • ncep_post.fd
  • OTLIFT.f
  • Generated by - doxygen 1.8.6
    • + doxygen 1.8.5 diff --git a/ZENSUN_8f.html b/ZENSUN_8f.html index 34fba4476..8854679c8 100644 --- a/ZENSUN_8f.html +++ b/ZENSUN_8f.html @@ -3,7 +3,7 @@ - + UPP: ZENSUN.f File Reference @@ -30,7 +30,7 @@
    UPP -  V11.0.0 +  11.0.0
    @@ -38,7 +38,7 @@ - + @@ -196,12 +196,13 @@

    This data is characterized by 74 points.

    Parameters
    - - - - - - + + + + + + +
    [in]dayJulian day (positive scalar or vector), (spring equinox = 80), (summer solstice= 171), (fall equinox = 266), (winter solstice= 356).
    [in]timeUniversal Time in hours (scalar or vector).
    [in]latGeographic latitude of point on earth's surface (degrees).
    [in]lonGeographic longitude of point on earth's surface (degrees).
    [out]sun_zenith- solar zenith angle.
    [out]sun_azimuth- solar azimuth angle.
    [in]dayinteger Julian day (positive scalar or vector), (spring equinox = 80), (summer solstice= 171), (fall equinox = 266), (winter solstice= 356).
    [in]timereal Universal Time in hours (scalar or vector).
    [in]latreal Geographic latitude of point on earth's surface (degrees).
    [in]lonreal Geographic longitude of point on earth's surface (degrees).
    [in]pireal The mathematical constant pi.
    [out]sun_zenithreal Solar zenith angle.
    [out]sun_azimuthreal Solar azimuth angle.
    @@ -215,9 +216,9 @@

    Program history log:

    Author
    Paul Ricchiazzi Earth Space Research Group,UCSB
    Date
    1992-10-23
    -

    Definition at line 61 of file ZENSUN.f.

    +

    Definition at line 62 of file ZENSUN.f.

    -

    Referenced by INITPOST().

    +

    Referenced by initpost().

    @@ -229,7 +230,7 @@

    Program history log:

  • sorc
  • ncep_post.fd
  • ZENSUN.f
  • Generated by - doxygen 1.8.6
    • + doxygen 1.8.5 diff --git a/ZENSUN_8f_source.html b/ZENSUN_8f_source.html index 5d5192612..d9419690f 100644 --- a/ZENSUN_8f_source.html +++ b/ZENSUN_8f_source.html @@ -3,7 +3,7 @@ - + UPP: ZENSUN.f Source File @@ -30,7 +30,7 @@
    UPP -  V11.0.0 +  11.0.0
    @@ -38,7 +38,7 @@ - + @@ -108,174 +108,174 @@
    Go to the documentation of this file.
    1 
    5 
    -
    61  subroutine zensun(day,time,lat,lon,pi,sun_zenith,sun_azimuth)
    -
    62 
    -
    63 !
    -
    64  use kinds, only: r_kind,i_kind
    -
    65 
    -
    66  implicit none
    -
    67 
    -
    68  integer(i_kind),intent(in):: day
    -
    69  real(r_kind),intent(in) :: pi,time,lat,lon
    -
    70  integer(i_kind) di
    -
    71  real(r_kind) z,gaz
    -
    72  real(r_kind) ut,noon
    -
    73  real(r_kind) y(5),y2(5),x(2,5),tx(5,2),xtx(2,2),atx(5,2),det
    -
    74  real(r_kind) tt,eqtime,decang,latsun,lonsun
    -
    75  real(r_kind) nday(74),eqt(74),dec(74)
    -
    76  real(r_kind) beta(2), beta2(2), a(2,2)
    -
    77  real(r_kind) t0,t1,p0,p1,zz,xx,yy
    -
    78  real(r_kind),intent(out) :: sun_zenith,sun_azimuth
    -
    79 
    -
    80  data nday/1.0, 6.0, 11.0, 16.0, 21.0, 26.0, 31.0, 36.0, 41.0, 46.0,&
    -
    81  51.0, 56.0, 61.0, 66.0, 71.0, 76.0, 81.0, 86.0, 91.0, 96.0,&
    -
    82  101.0, 106.0, 111.0, 116.0, 121.0, 126.0, 131.0, 136.0, 141.0, 146.0,&
    -
    83  151.0, 156.0, 161.0, 166.0, 171.0, 176.0, 181.0, 186.0, 191.0, 196.0,&
    -
    84  201.0, 206.0, 211.0, 216.0, 221.0, 226.0, 231.0, 236.0, 241.0, 246.0,&
    -
    85  251.0, 256.0, 261.0, 266.0, 271.0, 276.0, 281.0, 286.0, 291.0, 296.0,&
    -
    86  301.0, 306.0, 311.0, 316.0, 321.0, 326.0, 331.0, 336.0, 341.0, 346.0,&
    -
    87  351.0, 356.0, 361.0, 366.0/
    -
    88 
    -
    89  data eqt/ -3.23, -5.49, -7.60, -9.48,-11.09,-12.39,-13.34,-13.95,-14.23,-14.19,&
    -
    90  -13.85,-13.22,-12.35,-11.26,-10.01, -8.64, -7.18, -5.67, -4.16, -2.69,&
    -
    91  -1.29, -0.02, 1.10, 2.05, 2.80, 3.33, 3.63, 3.68, 3.49, 3.09,&
    -
    92  2.48, 1.71, 0.79, -0.24, -1.33, -2.41, -3.45, -4.39, -5.20, -5.84,&
    -
    93  -6.28, -6.49, -6.44, -6.15, -5.60, -4.82, -3.81, -2.60, -1.19, 0.36,&
    -
    94  2.03, 3.76, 5.54, 7.31, 9.04, 10.69, 12.20, 13.53, 14.65, 15.52,&
    -
    95  16.12, 16.41, 16.36, 15.95, 15.19, 14.09, 12.67, 10.93, 8.93, 6.70,&
    -
    96  4.32, 1.86, -0.62, -3.23/
    -
    97 
    -
    98  data dec/ -23.06,-22.57,-21.91,-21.06,-20.05,-18.88,-17.57,-16.13,-14.57,-12.91,&
    -
    99  -11.16, -9.34, -7.46, -5.54, -3.59, -1.62, 0.36, 2.33, 4.28, 6.19,&
    -
    100  8.06, 9.88, 11.62, 13.29, 14.87, 16.34, 17.70, 18.94, 20.04, 21.00,&
    -
    101  21.81, 22.47, 22.95, 23.28, 23.43, 23.40, 23.21, 22.85, 22.32, 21.63,&
    -
    102  20.79, 19.80, 18.67, 17.42, 16.05, 14.57, 13.00, 11.33, 9.60, 7.80,&
    -
    103  5.95, 4.06, 2.13, 0.19, -1.75, -3.69, -5.62, -7.51, -9.36,-11.16,&
    -
    104  -12.88,-14.53,-16.07,-17.50,-18.81,-19.98,-20.99,-21.85,-22.52,-23.02,&
    -
    105  -23.33,-23.44,-23.35,-23.06/
    -
    106 
    -
    107 !
    -
    108 ! compute the subsolar coordinates
    -
    109 !
    -
    110 
    -
    111  tt= mod(real((int(day)+time/24.-1.)),365.25) +1. ! fractional day number
    -
    112  ! with 12am 1jan = 1.
    -
    113  do di = 1, 73
    -
    114  if ((tt >= nday(di)) .and. (tt <= nday(di+1))) exit
    -
    115  end do
    -
    116 
    -
    117 !============== Perform a least squares regression on doy**3 ==============
    -
    118 
    -
    119  x(1,:) = 1.0
    -
    120 
    -
    121  if ((di >= 3) .and. (di <= 72)) then
    -
    122  y(:) = eqt(di-2:di+2)
    -
    123  y2(:) = dec(di-2:di+2)
    -
    124 
    -
    125  x(2,:) = nday(di-2:di+2)**3
    -
    126  end if
    -
    127  if (di == 2) then
    -
    128  y(1) = eqt(73)
    -
    129  y(2:5) = eqt(di-1:di+2)
    -
    130  y2(1) = dec(73)
    -
    131  y2(2:5) = dec(di-1:di+2)
    -
    132 
    -
    133  x(2,1) = nday(73)**3
    -
    134  x(2,2:5) = (365.+nday(di-1:di+2))**3
    -
    135  end if
    -
    136  if (di == 1) then
    -
    137  y(1:2) = eqt(72:73)
    -
    138  y(3:5) = eqt(di:di+2)
    -
    139  y2(1:2) = dec(72:73)
    -
    140  y2(3:5) = dec(di:di+2)
    -
    141 
    -
    142  x(2,1:2) = nday(72:73)**3
    -
    143  x(2,3:5) = (365.+nday(di:di+2))**3
    -
    144  end if
    -
    145  if (di == 73) then
    -
    146  y(1:4) = eqt(di-2:di+1)
    -
    147  y(5) = eqt(2)
    -
    148  y2(1:4) = dec(di-2:di+1)
    -
    149  y2(5) = dec(2)
    -
    150 
    -
    151  x(2,1:4) = nday(di-2:di+1)**3
    -
    152  x(2,5) = (365.+nday(2))**3
    -
    153  end if
    -
    154  if (di == 74) then
    -
    155  y(1:3) = eqt(di-2:di)
    -
    156  y(4:5) = eqt(2:3)
    -
    157  y2(1:3) = dec(di-2:di)
    -
    158  y2(4:5) = dec(2:3)
    -
    159 
    -
    160  x(2,1:3) = nday(di-2:di)**3
    -
    161  x(2,4:5) = (365.+nday(2:3))**3
    -
    162  end if
    -
    163 
    -
    164 ! Tx = transpose(x)
    -
    165  tx(1:5,1)=x(1,1:5)
    -
    166  tx(1:5,2)=x(2,1:5)
    -
    167 ! xTx = MATMUL(x,Tx)
    -
    168  xtx(1,1)=x(1,1)*tx(1,1)+x(1,2)*tx(2,1)+x(1,3)*tx(3,1)+x(1,4)*tx(4,1)+x(1,5)*tx(5,1)
    -
    169  xtx(1,2)=x(1,1)*tx(1,2)+x(1,2)*tx(2,2)+x(1,3)*tx(3,2)+x(1,4)*tx(4,2)+x(1,5)*tx(5,2)
    -
    170  xtx(2,1)=x(2,1)*tx(1,1)+x(2,2)*tx(2,1)+x(2,3)*tx(3,1)+x(2,4)*tx(4,1)+x(2,5)*tx(5,1)
    -
    171  xtx(2,2)=x(2,1)*tx(1,2)+x(2,2)*tx(2,2)+x(2,3)*tx(3,2)+x(2,4)*tx(4,2)+x(2,5)*tx(5,2)
    -
    172 
    -
    173  det = xtx(1,1)*xtx(2,2) - xtx(1,2)*xtx(2,1)
    -
    174  a(1,1) = xtx(2,2)/det
    -
    175  a(1,2) = -xtx(1,2)/det
    -
    176  a(2,1) = -xtx(2,1)/det
    -
    177  a(2,2) = xtx(1,1)/det
    -
    178 
    -
    179 ! aTx = MATMUL(Tx,a)
    -
    180  atx(1,1)=tx(1,1)*a(1,1)+tx(1,2)*a(2,1)
    -
    181  atx(2,1)=tx(2,1)*a(1,1)+tx(2,2)*a(2,1)
    -
    182  atx(3,1)=tx(3,1)*a(1,1)+tx(3,2)*a(2,1)
    -
    183  atx(4,1)=tx(4,1)*a(1,1)+tx(4,2)*a(2,1)
    -
    184  atx(5,1)=tx(5,1)*a(1,1)+tx(5,2)*a(2,1)
    -
    185  atx(1,2)=tx(1,1)*a(1,2)+tx(1,2)*a(2,2)
    -
    186  atx(2,2)=tx(2,1)*a(1,2)+tx(2,2)*a(2,2)
    -
    187  atx(3,2)=tx(3,1)*a(1,2)+tx(3,2)*a(2,2)
    -
    188  atx(4,2)=tx(4,1)*a(1,2)+tx(4,2)*a(2,2)
    -
    189  atx(5,2)=tx(5,1)*a(1,2)+tx(5,2)*a(2,2)
    -
    190 
    -
    191 ! beta = MATMUL(y,aTx)
    -
    192  beta(1) = y(1)*atx(1,1)+y(2)*atx(2,1)+y(3)*atx(3,1)+y(4)*atx(4,1)+y(5)*atx(5,1)
    -
    193  beta(2) = y(1)*atx(1,2)+y(2)*atx(2,2)+y(3)*atx(3,2)+y(4)*atx(4,2)+y(5)*atx(5,2)
    -
    194 
    -
    195 ! beta2 = MATMUL(y2,aTx)
    -
    196  beta2(1) = y2(1)*atx(1,1)+y2(2)*atx(2,1)+y2(3)*atx(3,1)+y2(4)*atx(4,1)+y2(5)*atx(5,1)
    -
    197  beta2(2) = y2(1)*atx(1,2)+y2(2)*atx(2,2)+y2(3)*atx(3,2)+y2(4)*atx(4,2)+y2(5)*atx(5,2)
    -
    198 
    -
    199 !============== finished least squares regression on doy**3 ==============
    -
    200 
    -
    201  if ((di < 3) .or. (di > 72)) tt = tt + 365.
    -
    202 
    -
    203  eqtime=(beta(1) + beta(2)*tt**3)/60.
    -
    204  decang=beta2(1) + beta2(2)*tt**3
    -
    205  latsun=decang
    -
    206 
    -
    207  ut=time
    -
    208  noon=12.-lon/15. ! universal time of noon
    -
    209 
    -
    210  lonsun=-15.*(ut-12.+eqtime)
    -
    211 
    -
    212  t0=(90.-lat)*pi/180.
    -
    213  t1=(90.-latsun)*pi/180.
    -
    214 
    -
    215  p0=lon*pi/180.
    -
    216  p1=lonsun*pi/180.
    -
    217 
    -
    218  zz=cos(t0)*cos(t1)+sin(t0)*sin(t1)*cos(p1-p0)
    -
    219 ! zz2=sin(t0)*sin(t1)+cos(t0)*cos(t1)*cos(p1-p0)
    -
    220  xx=sin(t1)*sin(p1-p0)
    -
    221  yy=sin(t0)*cos(t1)-cos(t0)*sin(t1)*cos(p1-p0)
    -
    222 
    -
    223  sun_zenith=90-acos(zz)/(pi/180)
    -
    224  sun_azimuth=atan2(xx,yy)/(pi/180)
    -
    225 
    -
    226  return
    -
    227 end subroutine zensun
    -
    zensun
    subroutine zensun(day, time, lat, lon, pi, sun_zenith, sun_azimuth)
    This subroutine computes solar position information as a function of geographic coordinates, date and time.
    Definition: ZENSUN.f:61
    +
    62  subroutine zensun(day,time,lat,lon,pi,sun_zenith,sun_azimuth)
    +
    63 
    +
    64 !
    +
    65  use kinds, only: r_kind,i_kind
    +
    66 
    +
    67  implicit none
    +
    68 
    +
    69  integer(i_kind),intent(in):: day
    +
    70  real(r_kind),intent(in) :: pi,time,lat,lon
    +
    71  integer(i_kind) di
    +
    72  real(r_kind) z,gaz
    +
    73  real(r_kind) ut,noon
    +
    74  real(r_kind) y(5),y2(5),x(2,5),tx(5,2),xtx(2,2),atx(5,2),det
    +
    75  real(r_kind) tt,eqtime,decang,latsun,lonsun
    +
    76  real(r_kind) nday(74),eqt(74),dec(74)
    +
    77  real(r_kind) beta(2), beta2(2), a(2,2)
    +
    78  real(r_kind) t0,t1,p0,p1,zz,xx,yy
    +
    79  real(r_kind),intent(out) :: sun_zenith,sun_azimuth
    +
    80 
    +
    81  data nday/1.0, 6.0, 11.0, 16.0, 21.0, 26.0, 31.0, 36.0, 41.0, 46.0,&
    +
    82  51.0, 56.0, 61.0, 66.0, 71.0, 76.0, 81.0, 86.0, 91.0, 96.0,&
    +
    83  101.0, 106.0, 111.0, 116.0, 121.0, 126.0, 131.0, 136.0, 141.0, 146.0,&
    +
    84  151.0, 156.0, 161.0, 166.0, 171.0, 176.0, 181.0, 186.0, 191.0, 196.0,&
    +
    85  201.0, 206.0, 211.0, 216.0, 221.0, 226.0, 231.0, 236.0, 241.0, 246.0,&
    +
    86  251.0, 256.0, 261.0, 266.0, 271.0, 276.0, 281.0, 286.0, 291.0, 296.0,&
    +
    87  301.0, 306.0, 311.0, 316.0, 321.0, 326.0, 331.0, 336.0, 341.0, 346.0,&
    +
    88  351.0, 356.0, 361.0, 366.0/
    +
    89 
    +
    90  data eqt/ -3.23, -5.49, -7.60, -9.48,-11.09,-12.39,-13.34,-13.95,-14.23,-14.19,&
    +
    91  -13.85,-13.22,-12.35,-11.26,-10.01, -8.64, -7.18, -5.67, -4.16, -2.69,&
    +
    92  -1.29, -0.02, 1.10, 2.05, 2.80, 3.33, 3.63, 3.68, 3.49, 3.09,&
    +
    93  2.48, 1.71, 0.79, -0.24, -1.33, -2.41, -3.45, -4.39, -5.20, -5.84,&
    +
    94  -6.28, -6.49, -6.44, -6.15, -5.60, -4.82, -3.81, -2.60, -1.19, 0.36,&
    +
    95  2.03, 3.76, 5.54, 7.31, 9.04, 10.69, 12.20, 13.53, 14.65, 15.52,&
    +
    96  16.12, 16.41, 16.36, 15.95, 15.19, 14.09, 12.67, 10.93, 8.93, 6.70,&
    +
    97  4.32, 1.86, -0.62, -3.23/
    +
    98 
    +
    99  data dec/ -23.06,-22.57,-21.91,-21.06,-20.05,-18.88,-17.57,-16.13,-14.57,-12.91,&
    +
    100  -11.16, -9.34, -7.46, -5.54, -3.59, -1.62, 0.36, 2.33, 4.28, 6.19,&
    +
    101  8.06, 9.88, 11.62, 13.29, 14.87, 16.34, 17.70, 18.94, 20.04, 21.00,&
    +
    102  21.81, 22.47, 22.95, 23.28, 23.43, 23.40, 23.21, 22.85, 22.32, 21.63,&
    +
    103  20.79, 19.80, 18.67, 17.42, 16.05, 14.57, 13.00, 11.33, 9.60, 7.80,&
    +
    104  5.95, 4.06, 2.13, 0.19, -1.75, -3.69, -5.62, -7.51, -9.36,-11.16,&
    +
    105  -12.88,-14.53,-16.07,-17.50,-18.81,-19.98,-20.99,-21.85,-22.52,-23.02,&
    +
    106  -23.33,-23.44,-23.35,-23.06/
    +
    107 
    +
    108 !
    +
    109 ! compute the subsolar coordinates
    +
    110 !
    +
    111 
    +
    112  tt= mod(real((int(day)+time/24.-1.)),365.25) +1. ! fractional day number
    +
    113  ! with 12am 1jan = 1.
    +
    114  do di = 1, 73
    +
    115  if ((tt >= nday(di)) .and. (tt <= nday(di+1))) exit
    +
    116  end do
    +
    117 
    +
    118 !============== Perform a least squares regression on doy**3 ==============
    +
    119 
    +
    120  x(1,:) = 1.0
    +
    121 
    +
    122  if ((di >= 3) .and. (di <= 72)) then
    +
    123  y(:) = eqt(di-2:di+2)
    +
    124  y2(:) = dec(di-2:di+2)
    +
    125 
    +
    126  x(2,:) = nday(di-2:di+2)**3
    +
    127  end if
    +
    128  if (di == 2) then
    +
    129  y(1) = eqt(73)
    +
    130  y(2:5) = eqt(di-1:di+2)
    +
    131  y2(1) = dec(73)
    +
    132  y2(2:5) = dec(di-1:di+2)
    +
    133 
    +
    134  x(2,1) = nday(73)**3
    +
    135  x(2,2:5) = (365.+nday(di-1:di+2))**3
    +
    136  end if
    +
    137  if (di == 1) then
    +
    138  y(1:2) = eqt(72:73)
    +
    139  y(3:5) = eqt(di:di+2)
    +
    140  y2(1:2) = dec(72:73)
    +
    141  y2(3:5) = dec(di:di+2)
    +
    142 
    +
    143  x(2,1:2) = nday(72:73)**3
    +
    144  x(2,3:5) = (365.+nday(di:di+2))**3
    +
    145  end if
    +
    146  if (di == 73) then
    +
    147  y(1:4) = eqt(di-2:di+1)
    +
    148  y(5) = eqt(2)
    +
    149  y2(1:4) = dec(di-2:di+1)
    +
    150  y2(5) = dec(2)
    +
    151 
    +
    152  x(2,1:4) = nday(di-2:di+1)**3
    +
    153  x(2,5) = (365.+nday(2))**3
    +
    154  end if
    +
    155  if (di == 74) then
    +
    156  y(1:3) = eqt(di-2:di)
    +
    157  y(4:5) = eqt(2:3)
    +
    158  y2(1:3) = dec(di-2:di)
    +
    159  y2(4:5) = dec(2:3)
    +
    160 
    +
    161  x(2,1:3) = nday(di-2:di)**3
    +
    162  x(2,4:5) = (365.+nday(2:3))**3
    +
    163  end if
    +
    164 
    +
    165 ! Tx = transpose(x)
    +
    166  tx(1:5,1)=x(1,1:5)
    +
    167  tx(1:5,2)=x(2,1:5)
    +
    168 ! xTx = MATMUL(x,Tx)
    +
    169  xtx(1,1)=x(1,1)*tx(1,1)+x(1,2)*tx(2,1)+x(1,3)*tx(3,1)+x(1,4)*tx(4,1)+x(1,5)*tx(5,1)
    +
    170  xtx(1,2)=x(1,1)*tx(1,2)+x(1,2)*tx(2,2)+x(1,3)*tx(3,2)+x(1,4)*tx(4,2)+x(1,5)*tx(5,2)
    +
    171  xtx(2,1)=x(2,1)*tx(1,1)+x(2,2)*tx(2,1)+x(2,3)*tx(3,1)+x(2,4)*tx(4,1)+x(2,5)*tx(5,1)
    +
    172  xtx(2,2)=x(2,1)*tx(1,2)+x(2,2)*tx(2,2)+x(2,3)*tx(3,2)+x(2,4)*tx(4,2)+x(2,5)*tx(5,2)
    +
    173 
    +
    174  det = xtx(1,1)*xtx(2,2) - xtx(1,2)*xtx(2,1)
    +
    175  a(1,1) = xtx(2,2)/det
    +
    176  a(1,2) = -xtx(1,2)/det
    +
    177  a(2,1) = -xtx(2,1)/det
    +
    178  a(2,2) = xtx(1,1)/det
    +
    179 
    +
    180 ! aTx = MATMUL(Tx,a)
    +
    181  atx(1,1)=tx(1,1)*a(1,1)+tx(1,2)*a(2,1)
    +
    182  atx(2,1)=tx(2,1)*a(1,1)+tx(2,2)*a(2,1)
    +
    183  atx(3,1)=tx(3,1)*a(1,1)+tx(3,2)*a(2,1)
    +
    184  atx(4,1)=tx(4,1)*a(1,1)+tx(4,2)*a(2,1)
    +
    185  atx(5,1)=tx(5,1)*a(1,1)+tx(5,2)*a(2,1)
    +
    186  atx(1,2)=tx(1,1)*a(1,2)+tx(1,2)*a(2,2)
    +
    187  atx(2,2)=tx(2,1)*a(1,2)+tx(2,2)*a(2,2)
    +
    188  atx(3,2)=tx(3,1)*a(1,2)+tx(3,2)*a(2,2)
    +
    189  atx(4,2)=tx(4,1)*a(1,2)+tx(4,2)*a(2,2)
    +
    190  atx(5,2)=tx(5,1)*a(1,2)+tx(5,2)*a(2,2)
    +
    191 
    +
    192 ! beta = MATMUL(y,aTx)
    +
    193  beta(1) = y(1)*atx(1,1)+y(2)*atx(2,1)+y(3)*atx(3,1)+y(4)*atx(4,1)+y(5)*atx(5,1)
    +
    194  beta(2) = y(1)*atx(1,2)+y(2)*atx(2,2)+y(3)*atx(3,2)+y(4)*atx(4,2)+y(5)*atx(5,2)
    +
    195 
    +
    196 ! beta2 = MATMUL(y2,aTx)
    +
    197  beta2(1) = y2(1)*atx(1,1)+y2(2)*atx(2,1)+y2(3)*atx(3,1)+y2(4)*atx(4,1)+y2(5)*atx(5,1)
    +
    198  beta2(2) = y2(1)*atx(1,2)+y2(2)*atx(2,2)+y2(3)*atx(3,2)+y2(4)*atx(4,2)+y2(5)*atx(5,2)
    +
    199 
    +
    200 !============== finished least squares regression on doy**3 ==============
    +
    201 
    +
    202  if ((di < 3) .or. (di > 72)) tt = tt + 365.
    +
    203 
    +
    204  eqtime=(beta(1) + beta(2)*tt**3)/60.
    +
    205  decang=beta2(1) + beta2(2)*tt**3
    +
    206  latsun=decang
    +
    207 
    +
    208  ut=time
    +
    209  noon=12.-lon/15. ! universal time of noon
    +
    210 
    +
    211  lonsun=-15.*(ut-12.+eqtime)
    +
    212 
    +
    213  t0=(90.-lat)*pi/180.
    +
    214  t1=(90.-latsun)*pi/180.
    +
    215 
    +
    216  p0=lon*pi/180.
    +
    217  p1=lonsun*pi/180.
    +
    218 
    +
    219  zz=cos(t0)*cos(t1)+sin(t0)*sin(t1)*cos(p1-p0)
    +
    220 ! zz2=sin(t0)*sin(t1)+cos(t0)*cos(t1)*cos(p1-p0)
    +
    221  xx=sin(t1)*sin(p1-p0)
    +
    222  yy=sin(t0)*cos(t1)-cos(t0)*sin(t1)*cos(p1-p0)
    +
    223 
    +
    224  sun_zenith=90-acos(zz)/(pi/180)
    +
    225  sun_azimuth=atan2(xx,yy)/(pi/180)
    +
    226 
    +
    227  return
    +
    228 end subroutine zensun
    +
    zensun
    subroutine zensun(day, time, lat, lon, pi, sun_zenith, sun_azimuth)
    This subroutine computes solar position information as a function of geographic coordinates, date and time.
    Definition: ZENSUN.f:62
    @@ -285,7 +285,7 @@
  • sorc
  • ncep_post.fd
  • ZENSUN.f
  • Generated by - doxygen 1.8.6
    • + doxygen 1.8.5