pkg/fizhi/getpwhere.F

a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
45f4a0300f Jean*0002 #undef CHECK_GETPWHERE
                0003
28e90e66e5 Andr*0004        subroutine getpwhere(myThid,numpress,pressures,levpressures)
45f4a0300f Jean*0005 C***********************************************************************
                0006 C subroutine getpwhere
                0007 C
                0008 C Purpose: Approximate (!) the level at which the mid-level pressure
                0009 C          is less than (ie, above in the atmosphere) a given value.
                0010 C
                0011 C Algorithm: Assume surface pressure of 1000 mb, and the pressure
bfa80ec297 Andr*0012 c          thicknesses set in make_phys_grid, with drF thicknesses above
45f4a0300f Jean*0013 C
                0014 C Need:    Information about the dynamics grid vertical spacing
                0015 C
                0016 C Input:   myThid       - process(or) number
                0017 C          numpres      - Number of pressures to process
                0018 C          pressures    - Pressure values to find levels for
                0019 C
                0020 C Output:  levpressures - Array of levels at which pressures are found
                0021 C                         These pressure levels correspond to the fizhi
                0022 C                         levels, and assume that the levels are counted
                0023 C                         from top to bottom (CRITICAL!)
                0024 C
                0025 C NOTE: The new physics levels specified here MUST correspond to the
                0026 C       physics levels specified in make_phys_grid from gridalt package.
                0027 C***********************************************************************
bfa80ec297 Andr*0028        implicit none
                0029 c
                0030 #include "SIZE.h"
                0031 #include "fizhi_SIZE.h"
                0032 #include "GRID.h"
0889f02121 Jean*0033 #ifdef CHECK_GETPWHERE
                0034 #include "EEPARAMS.h"
                0035 #include "PARAMS.h"
                0036 #endif /* CHECK_GETPWHERE */
bfa80ec297 Andr*0037
                0038        integer myThid,numpress
                0039        _RL pressures(numpress)
                0040        integer levpressures(numpress)
                0041 c
0889f02121 Jean*0042 #ifdef CHECK_GETPWHERE
                0043        INTEGER ioUnit, upperBnd
                0044 c      CHARACTER*(MAX_LEN_MBUF) msgBuf
                0045 #endif /* CHECK_GETPWHERE */
bfa80ec297 Andr*0046        integer n,L,dynlev
45f4a0300f Jean*0047 C Code that MUST correspond to make_phys_grid in the gridalt package!
                0048 C Require 12 bottom levels (300 mb worth) for the physics,
3daafce20b Jean*0049 C    Counting from bottom to top, the dp(s) are:
c3bac45152 Andr*0050        integer ntry,ntry10,ntry40
                0051        parameter (ntry40=15)
                0052        parameter (ntry10=12)
                0053        _RL dptry(ntry40), dptry10(ntry10), dptry40(ntry40)
                0054        _RL dptry_pedge(ntry40+1)
45f4a0300f Jean*0055 #ifdef TRY_NEW_GETPWHERE
                0056        _RL rC_dyn(Nr), dpTop
                0057 #else
c3bac45152 Andr*0058        _RL rF_pmid(Nr),rF_edge(Nr+1)
45f4a0300f Jean*0059 #endif
0889f02121 Jean*0060        _RL pref(Nrphys)
45f4a0300f Jean*0061        integer tmplev
c3bac45152 Andr*0062        data dptry10 /3.00, 6.00,10.00,14.00,17.00,25.00,
                0063      .              25.00,25.00,25.00,50.00,50.00,50.00/
                0064        data dptry40 /3.00, 6.00,10.00,14.00,17.00,25.00,
                0065      .              25.00,25.00,25.00,25.00,25.00,25.00,
                0066      .              25.00,25.00,25.00/
                0067
                0068        if( (Nr.eq.10).or.(Nr.eq.20) ) then
                0069         ntry = ntry10
                0070         do L = 1,ntry
                0071          dptry(L) = dptry10(L)
                0072         enddo
c909a9788d Andr*0073        elseif ((Nr.eq.40).or.(Nr.eq.70)) then
c3bac45152 Andr*0074         ntry = ntry40
                0075         do L = 1,ntry
                0076          dptry(L) = dptry40(L)
                0077         enddo
                0078        else
                0079         print *,' Dont know how to set levels for given pressures '
                0080         stop
                0081        endif
                0082
45f4a0300f Jean*0083 C define the mid pressure for the levels that are specified - bottom 300 mb.
63f7517fd8 Jean*0084 #ifdef TRY_NEW_GETPWHERE
                0085        dptry_pedge(1) = rF(1)*1. _d -2
                0086 #else  /* TRY_NEW_GETPWHERE */
c3bac45152 Andr*0087        dptry_pedge(1) = 1000.
63f7517fd8 Jean*0088 #endif /* TRY_NEW_GETPWHERE */
45f4a0300f Jean*0089        do L = 1,ntry
                0090         dptry_pedge(L+1) = dptry_pedge(L) - dptry(L)
c3bac45152 Andr*0091        enddo
                0092        do L = 1,ntry
                0093         pref(L) = (dptry_pedge(L) + dptry_pedge(L+1))/2.
                0094        enddo
45f4a0300f Jean*0095 #ifdef CHECK_GETPWHERE
                0096        ioUnit = errorMessageUnit
                0097        WRITE(ioUnit,'(A)') '===== GETPWHERE: CHECK start ====='
                0098        WRITE(ioUnit,'(4(A,I6),A)') '  Nr =', Nr,
                0099      &                        ' , Nrphys=', Nrphys,
                0100      &                        ' , ntry  =', ntry,
                0101      &                        ' , pref(1:ntry):'
                0102        WRITE(ioUnit,'(10F10.4)') (pref(L),L=1,ntry)
                0103 #endif /* CHECK_GETPWHERE */
                0104
                0105 #ifdef TRY_NEW_GETPWHERE
                0106 C top levels DP of 1 mb (or 0.01 mb for strat version)
                0107        dpTop = 1. _d 0
                0108        IF (Nr.EQ.70) dpTop = 1. _d -2
                0109        DO L = ntry+1,Nrphys
                0110         pref(L) = (Nrphys-L+0.5)*dpTop
                0111        ENDDO
                0112 C define the rest of the mid pressures from the dynamics levels
                0113        DO L = 1,Nr
                0114         rC_dyn(L) = rC(L)*1. _d -2
                0115        ENDDO
bfa80ec297 Andr*0116
45f4a0300f Jean*0117        dynlev = 0
                0118        DO L = 1,Nr
                0119          IF ( rC_dyn(L).GE.dptry_pedge(ntry+1) ) dynlev = L
                0120        ENDDO
                0121        DO L = ntry+1,MIN(Nrphys,ntry+Nr-dynlev)
                0122          pref(L) = rC_dyn(dynlev+L-ntry)
                0123        ENDDO
                0124 #else  /* TRY_NEW_GETPWHERE */
                0125 C define the rest of the mid pressures from the dynamics levels
c3bac45152 Andr*0126        rF_edge(1) = 1000.
                0127        do L = 2,Nr+1
                0128         rF_edge(L) = rF_edge(L-1) - (drF(L-1)/100.)
                0129        enddo
                0130        do L = 1,Nr
                0131         rF_pmid(L) = (rF_edge(L) + rF_edge(L+1))/2.
                0132        enddo
bfa80ec297 Andr*0133
c3bac45152 Andr*0134        dynlev = 0
                0135        do L = 1,Nr
                0136         if(rF_pmid(L).ge.pref(ntry)) dynlev = L
                0137        enddo
0889f02121 Jean*0138 #ifdef CHECK_GETPWHERE
                0139        IF ( rF_pmid(dynlev).ge.pref(ntry)-25. ) THEN
                0140          upperBnd = ntry+Nr-dynlev
                0141        ELSE
                0142          upperBnd = ntry+Nr-dynlev-1
                0143        ENDIF
45f4a0300f Jean*0144        WRITE(ioUnit,'(1(A,I5),A)') ' Up-Bnd=', upperBnd,
                0145      &                           ' , rF_pmid:'
0889f02121 Jean*0146        WRITE(ioUnit,'(10F10.4)') rF_pmid
45f4a0300f Jean*0147        IF ( upperBnd.LT.Nrphys-1 ) THEN
                0148          WRITE(ioUnit,'(A)')
                0149      &   'ERROR: exeeding "rF_pmid" array bounds => pref ill defined'
                0150 c        STOP 'ABNORMAL END: S/R GETPWHERE'
                0151        ENDIF
0889f02121 Jean*0152 #endif /* CHECK_GETPWHERE */
c3bac45152 Andr*0153        if(rF_pmid(dynlev).ge.pref(ntry)-25.) then
                0154         do L = ntry+1,Nrphys-1
                0155          pref(L) = rF_pmid(dynlev+L-ntry)
                0156         enddo
                0157        else
                0158         pref(ntry) = rF_pmid(dynlev)
                0159         do L = ntry+1,Nrphys-1
                0160          pref(L) = rF_pmid(dynlev+L-ntry+1)
                0161         enddo
                0162        endif
45f4a0300f Jean*0163 C Add top level DP of 1 mb - p mid is at 0.5 mb (or 0.05 mb for strat version)
c3bac45152 Andr*0164        pref(Nrphys) = 0.5
c909a9788d Andr*0165        if(Nr.eq.70)pref(Nrphys) = 0.005
45f4a0300f Jean*0166 #endif /* TRY_NEW_GETPWHERE */
                0167
                0168        DO n = 1,numpress
                0169         DO L = 1,Nrphys
                0170          IF ( pref(L).GE.pressures(n) ) tmplev = L
                0171         ENDDO
                0172
                0173 C and now flip the level numbers for the top down counting in fizhi
                0174         levpressures(n) = Nrphys+1-tmplev
                0175        ENDDO
bfa80ec297 Andr*0176
869d786604 Jean*0177 #ifdef CHECK_GETPWHERE
45f4a0300f Jean*0178        WRITE(ioUnit,'(3(A,I5),A)') ' dynlev=', dynlev,
                0179      &                           ' , numpress=', numpress,
                0180      &                           ' , pressures:'
869d786604 Jean*0181        WRITE(ioUnit,'(10F10.4)') pressures
45f4a0300f Jean*0182        WRITE(ioUnit,'(A)') 'pref(ntry:Nrphys):'
                0183        WRITE(ioUnit,'(10F10.4)') (pref(L),L=ntry,Nrphys)
869d786604 Jean*0184        WRITE(ioUnit,'(A)') 'levpressures:'
                0185        WRITE(ioUnit,'(20I5)') levpressures
                0186        WRITE(ioUnit,'(A)') '===== GETPWHERE: CHECK end   ====='
                0187 #endif /* CHECK_GETPWHERE */
                0188
45f4a0300f Jean*0189        RETURN
                0190        END