Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:40:09 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
8598abfc35 Andr*0002       SUBROUTINE QCHECK ( idim1,idim2,jdim1,jdim2,ldim,Nsx,Nsy,im1,im2,
                0003      .     jm1,jm2,bi,bj,dp,qz)
c4003ccfce Andr*0004 C***********************************************************************
337c673391 Jean*0005 C  Purpose
                0006 C     Check Specific Humidity Field for Negative values
c4003ccfce Andr*0007 C
                0008 C  Argument Description
e5b3b476cb Andr*0009 C     IDIM1 .... Start Zonal Dimension
                0010 C     IDIM2 .... End Zonal Dimension
                0011 C     JDIM1 .... Start Meridional Dimension
                0012 C     JDIM2 .... End Meridional Dimension
                0013 C     IM1 ...... Start Zonal Span
                0014 C     IM2 ...... End Zonal Span
                0015 C     JM1 ...... Start Meridional Span
                0016 C     JM2 ...... End Meridional Span
d10f990686 Andr*0017 C     LDIM ..... Vertical   Dimension
e5b3b476cb Andr*0018 C     DP ....... Delta Pressure
                0019 C     QZ ........Specific Humidity (g/g)
c4003ccfce Andr*0020 C***********************************************************************
                0021       implicit none
                0022 
e5b3b476cb Andr*0023       integer idim1,idim2,jdim1,jdim2,Ldim,im1,im2,jm1,jm2
8598abfc35 Andr*0024       integer Nsx,Nsy,bi,bj
                0025       _RL qz(idim1:idim2,jdim1:jdim2,ldim,Nsx,Nsy)
                0026       _RL dp(idim1:idim2,jdim1:jdim2,ldim,Nsx,Nsy)
c4003ccfce Andr*0027 
                0028       integer i,j,L,LM1
8598abfc35 Andr*0029       _RL ddsig
c4003ccfce Andr*0030 
                0031 c Fill Negative Specific Humidities
                0032 c ---------------------------------
e5b3b476cb Andr*0033       DO L=2,Ldim
337c673391 Jean*0034       LM1 = L-1
e5b3b476cb Andr*0035       do j=jm1,jm2
                0036       do i=im1,im2
337c673391 Jean*0037        ddsig = dp(i,j,LM1,bi,bj)/dp(i,j,L,bi,bj)
                0038        if( qz(i,j,LM1,bi,bj).lt.0.0  _d 0) then
8598abfc35 Andr*0039         qz(i,j,L,bi,bj  ) = qz(i,j,L,bi,bj) + qz(i,j,LM1,bi,bj)*ddsig
                0040         qz(i,j,LM1,bi,bj) = 0.0 _d 0
3d38dd8e43 Andr*0041        endif
                0042       enddo
                0043       enddo
                0044       enddo
c4003ccfce Andr*0045 
e5b3b476cb Andr*0046       do j=jm1,jm2
                0047       do i=im1,im2
8598abfc35 Andr*0048        if(qz(i,j,Ldim,bi,bj).lt.0.0  _d 0)qz(i,j,Ldim,bi,bj) = 0.0 _d 0
3d38dd8e43 Andr*0049       enddo
                0050       enddo
c4003ccfce Andr*0051 
337c673391 Jean*0052       return
                0053       end
c4003ccfce Andr*0054 
9e712e1a7f Andr*0055       subroutine tracer_fill ( pq,im,jm,lm,dlam,dphi,dp)
c4003ccfce Andr*0056 C***********************************************************************
337c673391 Jean*0057 C  PURPOSE
c4003ccfce Andr*0058 C     Fill negative tracer values using local borrowing
337c673391 Jean*0059 C
                0060 C  INPUT
c4003ccfce Andr*0061 C     pq ..... Mass-weighted (PI) Tracer
                0062 C     im ..... Zonal      Dimension
                0063 C     jm ..... Meridional Dimension
                0064 C     lm ..... Vertical   Dimension
                0065 C     dlam ... Zonal      Grid Increment
                0066 C     dphi ... Meridional Grid Increment
9e712e1a7f Andr*0067 C     dp ..... Vertical   Grid Increment
337c673391 Jean*0068 C
c4003ccfce Andr*0069 C  Note:
337c673391 Jean*0070 C     If no immediate surrounding value is large enough to fill negative
9e712e1a7f Andr*0071 C     value,
337c673391 Jean*0072 C     the sum of immediate surrounding positive values is tried.
c4003ccfce Andr*0073 C     If sum is not large enough, tracer is simply set to zero.
337c673391 Jean*0074 C
c4003ccfce Andr*0075 C***********************************************************************
                0076 
                0077       implicit none
                0078 
                0079 c Input Variables
                0080 c ---------------
                0081       integer im,jm,lm
a456aa407c Andr*0082       _RL    pq(im,jm,lm),dlam(im),dphi(jm),dp(im,jm,lm)
c4003ccfce Andr*0083 
                0084 c Local Variables
                0085 c ---------------
                0086       integer  i,j,l,im1,ip1,imax,m
a456aa407c Andr*0087       _RL     lam(im), phi(jm)
                0088       _RL     array(6)
                0089       _RL     pi,a,getcon,undef
                0090       _RL     qmax,qval,sum,fact
c4003ccfce Andr*0091 
a456aa407c Andr*0092       _RL        dxu(im,jm)
                0093       _RL        dxv(im,jm)
                0094       _RL        dxp(im,jm)
                0095       _RL        dyv(im,jm)
                0096       _RL        dyp(im,jm)
c4003ccfce Andr*0097 
a456aa407c Andr*0098       _RL d2p(im,jm)
c4003ccfce Andr*0099 
                0100 C *********************************************************
                0101 C ****                 Initialization                  ****
                0102 C *********************************************************
                0103 
                0104       pi = 4.0*atan(1.0)
                0105       a  = getcon('EARTH RADIUS')
                0106 
                0107 c Compute Longitudes
                0108 c ------------------
                0109       lam(1) = -pi
                0110       do i=2,im
                0111       lam(i) = lam(i-1) + dlam(i-1)
                0112       enddo
                0113 
                0114 c Compute Latitudes
                0115 c -----------------
                0116       phi(1) = -pi/2.
                0117       do j=2,jm-1
                0118       phi(j) = phi(j-1) + dphi(j-1)
                0119       enddo
                0120       phi(jm) =  pi/2.
                0121 
                0122 c Compute DXU and DYV
                0123 c -------------------
                0124       do j=2,jm-1
                0125       do i=1,im
                0126       dxu(i,j) = a*cos(phi(j))*dlam(i)
                0127       enddo
                0128       enddo
                0129 
                0130       do j=2,jm-2
                0131       do i=1,im
                0132       dyv(i,j) = a*dphi(j)
                0133       enddo
                0134       enddo
                0135       do i=1,im
                0136       dyv(i,1)    = a*(dphi(1)   +0.5*dphi(2)   )
                0137       dyv(i,jm-1) = a*(dphi(jm-1)+0.5*dphi(jm-2))
                0138       enddo
                0139 
                0140 c Compute DXP and DXV
                0141 c -------------------
                0142       do j=2,jm-1
                0143       im1 =  im
                0144       do i=1,im
                0145       dxp(i,j) = ( dxu(i,j)+dxu(im1,j) )*0.5
                0146       im1 = i
                0147       enddo
                0148       enddo
                0149 
                0150       do j=2,jm-2
                0151       do i=1,im
                0152       dxv(i,j) = ( dxp(i,j)+dxp(i,j+1) )*0.5
                0153       enddo
                0154       enddo
                0155 
                0156 c Compute DYP
                0157 c -----------
                0158       do j=3,jm-2
                0159       do i=1,im
                0160       dyp(i,j) = ( dyv(i,j)+dyv(i,j-1) )*0.5
                0161       enddo
                0162       enddo
                0163       do i=1,im
                0164       dyp(i,2)    = dyv(i,1)
                0165       dyp(i,jm-1) = dyv(i,jm-1)
                0166       enddo
                0167 
                0168 c Compute Area Factor D2P
                0169 c -----------------------
                0170       do j=3,jm-2
                0171       do i=1,im
                0172       d2p(i,j) = 0.5*( dxv(i,j)+dxv(i,j-1) )*dyp(i,j)
                0173       enddo
                0174       enddo
                0175       do i=1,im
                0176       d2p(i,2)    = dxv(i,2)   *dyp(i,2)
                0177       d2p(i,jm-1) = dxv(i,jm-2)*dyp(i,jm-1)
                0178       enddo
                0179 
                0180       undef = getcon('UNDEF')
                0181 
                0182 C *********************************************************
                0183 C ****             Fill Negative Values                ****
                0184 C *********************************************************
                0185 
                0186       do l=1,lm
                0187       do j=2,jm-1
                0188 
                0189       im1 = im-1
                0190       i   = im
                0191       do ip1=1,im
                0192 
                0193       if( pq(i,j,L).lt.0.0 ) then
                0194 
9e712e1a7f Andr*0195       qval     = pq(i  ,j,L)*d2p(i  ,j)*dp(i,j,L)
                0196       array(1) = pq(ip1,j,L)*d2p(ip1,j)*dp(i,j,L)
                0197       array(2) = pq(im1,j,L)*d2p(im1,j)*dp(i,j,L)
c4003ccfce Andr*0198 
                0199       if( j.eq.jm-1 ) then
                0200       array(3) = -undef
                0201       else
9e712e1a7f Andr*0202       array(3) = pq(i,j+1,L)*d2p(i,j+1)*dp(i,j,L)
c4003ccfce Andr*0203       endif
                0204       if( j.eq.2    ) then
                0205       array(4) = -undef
                0206       else
9e712e1a7f Andr*0207       array(4) = pq(i,j-1,L)*d2p(i,j-1)*dp(i,j,L)
c4003ccfce Andr*0208       endif
                0209       if( L.eq.1    ) then
                0210       array(5) = -undef
                0211       else
ed0b0d8f16 Andr*0212       array(5) = pq(i,j,L-1)*d2p(i,j)*dp(i,j,L)
c4003ccfce Andr*0213       endif
                0214       if( L.eq.lm   ) then
                0215       array(6) = -undef
                0216       else
ed0b0d8f16 Andr*0217       array(6) = pq(i,j,L+1)*d2p(i,j)*dp(i,j,L)
c4003ccfce Andr*0218       endif
                0219 
ed0b0d8f16 Andr*0220       call maxval1 (array,6,-qval,qmax,imax)
c4003ccfce Andr*0221 
                0222       if( imax.eq.0 ) then
                0223           sum = 0.0
                0224           do m=1,6
                0225           if( array(m).gt.0.0 ) sum = sum + array(m)
                0226           enddo
                0227              if( sum.gt.-qval ) then
                0228                fact = 1.0 + qval/sum
                0229                if( array(1).gt.0 ) pq(ip1,j,L) = pq(ip1,j,L) * fact
                0230                if( array(2).gt.0 ) pq(im1,j,L) = pq(im1,j,L) * fact
                0231                if( array(3).gt.0 ) pq(i,j+1,L) = pq(i,j+1,L) * fact
                0232                if( array(4).gt.0 ) pq(i,j-1,L) = pq(i,j-1,L) * fact
                0233                if( array(5).gt.0 ) pq(i,j,L-1) = pq(i,j,L-1) * fact
                0234                if( array(6).gt.0 ) pq(i,j,L+1) = pq(i,j,L+1) * fact
                0235                                    pq(i,j,L)   = 0.0
                0236              else
                0237                pq(i,j,L) = 0.0
                0238              endif
                0239       else
337c673391 Jean*0240           if( imax.eq.1 ) pq(ip1,j,L) = pq(ip1,j,L) +
9e712e1a7f Andr*0241      .                                 pq(i,j,L)*d2p(i,j)/d2p(ip1,j)
337c673391 Jean*0242           if( imax.eq.2 ) pq(im1,j,L) = pq(im1,j,L) +
9e712e1a7f Andr*0243      .                                 pq(i,j,L)*d2p(i,j)/d2p(im1,j)
337c673391 Jean*0244           if( imax.eq.3 ) pq(i,j+1,L) = pq(i,j+1,L) +
9e712e1a7f Andr*0245      .                                 pq(i,j,L)*d2p(i,j)/d2p(i,j+1)
337c673391 Jean*0246           if( imax.eq.4 ) pq(i,j-1,L) = pq(i,j-1,L) +
9e712e1a7f Andr*0247      .                                 pq(i,j,L)*d2p(i,j)/d2p(i,j-1)
337c673391 Jean*0248           if( imax.eq.5 ) pq(i,j,L-1) = pq(i,j,L-1) +
9e712e1a7f Andr*0249      .                                 pq(i,j,L)*dp(i,j,L) /dp(i,j,L-1)
337c673391 Jean*0250           if( imax.eq.6 ) pq(i,j,L+1) = pq(i,j,L+1) +
9e712e1a7f Andr*0251      .                                 pq(i,j,L)*dp(i,j,L) /dp(i,j,L+1)
c4003ccfce Andr*0252                           pq(i,j,L)   = 0.0
                0253       endif
                0254 
                0255       endif  ! End pq<0 Test
                0256 
                0257       im1 = i
                0258       i   = ip1
                0259       enddo
                0260       enddo
                0261       enddo
                0262 
                0263       return
                0264       end
                0265 
ed0b0d8f16 Andr*0266       subroutine maxval1 (q,im,qval,qmax,imax)
c4003ccfce Andr*0267 C***********************************************************************
337c673391 Jean*0268 C  PURPOSE
c4003ccfce Andr*0269 C     Find the location and value of the array element which is greater
                0270 C     than a prescribed value.
337c673391 Jean*0271 C
                0272 C  INPUT
                0273 C     q ...... Array Elements
                0274 C     im ..... Dimension of Array q
c4003ccfce Andr*0275 C     qval ... Prescribed Value
337c673391 Jean*0276 C
                0277 C  OUTPUT
                0278 C     qmax ... Largest Array element which is greater than qval
c4003ccfce Andr*0279 C     imax ... Location of Largest Array Element
337c673391 Jean*0280 C
c4003ccfce Andr*0281 C  Note:
337c673391 Jean*0282 C     If no array element is larger than qval, then imax = 0
                0283 C
c4003ccfce Andr*0284 C***********************************************************************
                0285       implicit none
                0286       integer  im, i, imax
a456aa407c Andr*0287       _RL   q(im), qmax, qval
c4003ccfce Andr*0288       qmax = qval
                0289       imax = 0
                0290       do i=1,im
                0291       if( q(i).gt.qmax ) then
                0292       qmax = q(i)
                0293       imax =   i
                0294       endif
                0295       enddo
                0296       return
                0297       end
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!