cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c This is the statistical convection parameterization c (Pasquero and Tziperman, JPO 2006). c c INPUT c NL : number of vertical layers c avgTMP : array of NL mean potential temperatures (deg C) c avgS : array of NL mean salinities (psu) c sigmaTMP: array of NL std deviations for temperatures (deg C) c sigmaS : array of NL std deviations for salinity (psu) c corr : array of NL TMP,S correlations ([-1,1]) c pres : array of NL pressures corresponding to each layer (dbar) c c OUTPUT c avgTMP : array of NL mean pot temp after convective adjustment c avgS : array of NL mean salinities after convective adjustment c sigmaTMP: array of NL std deviations for pot temperature c sigmaS : array of NL std deviations for salinity c corr : array of NL TMP,S corr after convective adjustment c cconv : array of NL values of probability of convection c before adjustment c c N.B.: Libraries for error function (derf) must be installed. c SUBROUTINE STAT_CONV_PARAM (NL,avgTMP,avgS,sigmaTMP + ,sigmaS,corr,pres,cconv) implicit none integer NL,i,j,irepeat,iupo,idowno,iup,idown,il,k integer iiter,icheck,N parameter(N=1000) real*8 avgTMP(N),avgS(N),sigmaTMP(N),sigmaS(N),corr(N) real*8 pres(N),conv,dz1,dz2,pr(0:N+1),sigma(N) real*8 sigma1,sigma2,avgt,avgb real*8 Pconv,Pnoconv,stot,den,expo,avg,argo real*8 q1,q2,q3,q4,dzinv,dz2inv,rh1,rh2,prn,Tnp real*8 avg1t,avg2t,avg1s,avg2s,cs1,cs2,sigma1t,sigma2t real*8 sigma1s,sigma2s,ang1t,ang2t,ang1s,ang2s,ang1tt,ang2tt real*8 ang1ss,ang2ss,ang1ts,ang2ts,ang1t2s,ang2t1s,ang12t,ang12s real*8 til1t,til2t,til1s,til2s,til1ss,til2ss,til1tt,til2tt real*8 til1ts,til2ts,tmixed,smixed,ttmixed,ssmixed,tsmixed real*8 cconv(N),alfa(N),beta(N),t0(N),s0(N) real*8 pi,soglia,cflag(N) parameter(pi=3.1415926535897932384626433,soglia=0.01) c soglia is the threshold in P_conv for convection to occur. common /leos/alfa,beta,t0,s0 c----------------------- c Initialize variables. c----------------------- do il=1,NL pr(il)=pres(il) if (sigmaTMP(il).eq.0.) sigmaTMP(il)=1.d-8 if (sigmaS(il).eq.0.) sigmaS(il)=1.d-8 cflag(il)=0. cconv(il)=0. enddo pr(0)=0. pr(NL+1)=2*pr(NL)-pr(NL-1) iupo=1 idowno=NL-1 c----------------------------------------------- c Linearize equation of state at each interface. c----------------------------------------------- do k=1,NL prn=(pr(k+1)+pr(k))*0.5 !******lower interface pressure call POTMP(0.d0,avgTMP(k),avgS(k),prn,Tnp) call RHO(prn,prn,(Tnp-.5),avgS(k),rh1) call RHO(prn,prn,(Tnp+.5),avgS(k),rh2) alfa(k)=-(rh2-rh1) call RHO(prn,prn,Tnp,(avgS(k)-0.2),rh1) call RHO(prn,prn,Tnp,(avgS(k)+0.2),rh2) beta(k)=(rh2-rh1)/(0.4) t0(k)=avgTMP(k) s0(k)=avgS(k) enddo c---------------------------------------------------------------------------------- c Instability check. c irepeat: loop over the vertical column, until everything is statistically stable c il : loop over the couples of nearby levels. c iiter : iterations of the procedure between two nearby levels, until statistical c stability is reached. c---------------------------------------------------------------------------------- do irepeat=1,100 !Repeat checks for new instabilities over the column. icheck=0 iup=iupo idown=idowno iupo=0 idowno=NL-1 do il=iup,idown ! Choose two nearby layer and checks for instability. dz1=(pr(il+1)-pr(il-1))*0.5 dz2=(pr(il+2)-pr(il))*0.5 dzinv=1./(dz1+dz2) dz2inv=1./(dz1+dz2)**2 c print *,dz1,dz2 iiter=0 10 continue iiter=iiter+1 call RHO_LIN(avgTMP(il),avgS(il),il,avgt) call RHO_LIN(avgTMP(il+1),avgS(il+1),il,avgb) avg1t=avgTMP(il) avg2t=avgTMP(il+1) avg1s=avgS(il) avg2s=avgS(il+1) sigma1t=sigmaTMP(il) sigma2t=sigmaTMP(il+1) sigma1s=sigmaS(il) sigma2s=sigmaS(il+1) cs1=corr(il) cs2=corr(il+1) c print *,"" c print *,'starting from mean T, S (top):',avg1t,avg1s c print *,'starting from mean T, S (bottom):',avg2t,avg2s c print *,'sigma T,S (top):',sigma1t,sigma1s c print *,'sigma T,S (bottom):',sigma2t,sigma2s c print *,'correl TS (top, bottom):',cs1,cs2 c write (10,*) avg1t,avg2t c write (10,*) avg1s,avg2s c write (10,*) sigma1t,sigma2t c write (10,*) sigma1s,sigma2s c write (10,*) cs1,cs2 sigma1=alfa(il)*alfa(il)*sigma1t**2+ + beta(il)**2*sigma1s**2-2.*alfa(il) + *beta(il)*sigma1t*sigma1s*cs1 sigma2=alfa(il)*alfa(il)*sigma2t**2+ + beta(il)**2*sigma2s**2-2.*alfa(il) + *beta(il)*sigma2t*sigma2s*cs2 stot=sigma1+sigma2 argo=(avgt-avgb)/dsqrt(2.*stot) Pconv=0.5*(1.d0+derf(argo)) !Probability of convection. c print *,'Pconv= ',Pconv if (iiter.eq.1.and.irepeat.eq.1) cconv(il)=Pconv c-------------- c Instability. c-------------- if (Pconv.gt.soglia) then icheck=1 iupo=max(iupo,il) idowno=min(idowno,il) expo=exp(-argo*argo) den=1./(dsqrt(2.*pi*stot)*Pconv) q1=-cs1*alfa(il)*sigma1t+beta(il)*sigma1s q2=-cs2*alfa(il)*sigma2t+beta(il)*sigma2s q3=-alfa(il)*sigma1t+cs1*beta(il)*sigma1s q4=-alfa(il)*sigma2t+cs2*beta(il)*sigma2s avg=(avgt-avgb)/stot ang1t=avg1t+sigma1t*den*expo*q3 ang1s=avg1s+sigma1s*den*expo*q1 ang2t=avg2t-sigma2t*den*expo*q4 ang2s=avg2s-sigma2s*den*expo*q2 ang1ts=(ang1s*avg1t+ang1t*avg1s-avg1s*avg1t)- + den*sigma1t*sigma1s*expo*avg*q3*q1+ + cs1*sigma1t*sigma1s ang2ts=(ang2s*avg2t+ang2t*avg2s-avg2s*avg2t)- + den*sigma2t*sigma2s*expo*avg*q4*q2+ + cs2*sigma2t*sigma2s ang1ss=avg1s**2+sigma1s**2+sigma1s*q1*den*expo* + (2.*avg1s-sigma1s*q1*avg) ang1tt=avg1t**2+sigma1t**2+sigma1t*q3*den*expo* + (2.*avg1t-sigma1t*q3*avg) ang2ss=avg2s**2+sigma2s**2-sigma2s*q2*den*expo* + (2.*avg2s+sigma2s*q2*avg) ang2tt=avg2t**2+sigma2t**2-sigma2t*q4*den*expo* + (2.*avg2t+sigma2t*q4*avg) ang1t2s=avg1t*ang2s+avg2s*ang1t-avg1t*avg2s+sigma1t* + sigma2s*den*expo*avg*q2*q3 ang2t1s=avg2t*ang1s+avg1s*ang2t-avg2t*avg1s+sigma2t* + sigma1s*den*expo*avg*q1*q4 ang12t=avg2t*ang1t+avg1t*ang2t-avg2t*avg1t+sigma2t* + sigma1t*den*expo*avg*q3*q4 ang12s=avg2s*ang1s+avg1s*ang2s-avg2s*avg1s+sigma2s* + sigma1s*den*expo*avg*q1*q2 Pnoconv=1.-Pconv if (Pnoconv.ne.0.) then den=1./(sqrt(2.*pi*stot)*Pnoconv) til1t=avg1t-sigma1t*den*expo*q3 til1s=avg1s-sigma1s*den*expo*q1 til2t=avg2t+sigma2t*den*expo*q4 til2s=avg2s+sigma2s*den*expo*q2 til1ts=(til1s*avg1t+til1t*avg1s-avg1s*avg1t)+ + den*sigma1t*sigma1s*expo*avg*q1*q3+ + cs1*sigma1t*sigma1s til2ts=(til2s*avg2t+til2t*avg2s-avg2s*avg2t)+ + den*sigma2t*sigma2s*expo*avg*q2*q4+ + cs2*sigma2t*sigma2s til1ss=avg1s**2+sigma1s**2-sigma1s*q1*den*expo* + (2.*avg1s-sigma1s*q1*avg) til1tt=avg1t**2+sigma1t**2-sigma1t*q3*den*expo* + (2.*avg1t-sigma1t*q3*avg) til2ss=avg2s**2+sigma2s**2+sigma2s*q2*den*expo* + (2.*avg2s+sigma2s*q2*avg) til2tt=avg2t**2+sigma2t**2+sigma2t*q4*den*expo* + (2.*avg2t+sigma2t*q4*avg) endif c print *,"Mean T,S, convective area:" c print *,' top=',real(ang1t),real(ang1s) c print *,' bottom=',real(ang2t),real(ang2s) c print *,"Mean T,S, non convective area:" c print *,' top=',real(til1t),real(til1s) c print *,' bottom=',real(til2t),real(til2s) c print *,"rms T,S, convective area:" c print *,' top=',sqrt(ang1tt-ang1t**2), c + sqrt(ang1ss-ang1s**2) c print *,' bottom=',sqrt(ang2tt-ang2t**2), c + sqrt(ang2ss-ang2s**2) c print *,"rms T,S, non convective area:" c print *,' top=',sqrt(til1tt-til1t**2), c + sqrt(til1ss-til1s**2) c print *,' bottom=',sqrt(til2tt-til2t**2), c + sqrt(til2ss-til2s**2) c print *,"Correlation TS, top and bottom, convective area" c print *,ang1ts,ang1t*ang1s c print *,(ang1ts-ang1t*ang1s)/(sqrt(ang1tt-ang1t**2)* c + sqrt(ang1ss-ang1s**2)),(ang2ts-ang2t*ang2s)/( c + sqrt(ang2tt-ang2t**2)*sqrt(ang2ss-ang2s**2)) c print *,"Correlation TS, top and bottom, nonconvective area" c print *,(til1ts-til1t*til1s)/(sqrt(til1tt-til1t**2)* c + sqrt(til1ss-til1s**2)),(til2ts-til2t*til2s)/ c + (sqrt(til2tt-til2t**2)*sqrt(til2ss-til2s**2)) c c Calculates mixed properties over convecting area c tmixed=(dz1*ang1t+dz2*ang2t)*dzinv smixed=(dz1*ang1s+dz2*ang2s)*dzinv ttmixed=(dz1*dz1*ang1tt+dz2*dz2*ang2tt+ + 2.*dz1*dz2*ang12t)*dz2inv ssmixed=(dz1*dz1*ang1ss+dz2*dz2*ang2ss+ + 2.*dz1*dz2*ang12s)*dz2inv tsmixed=(dz1*dz1*ang1ts+dz2*dz2*ang2ts+ + dz1*dz2*(ang1t2s+ang2t1s))*dz2inv c print *,'Mean T and S, mixed:',tmixed,smixed c print *,'Sigma T and S, mixed:', c + sqrt(ttmixed-tmixed**2), c + sqrt(ssmixed-smixed**2) c print *,'Correlation TS, mixed',(tsmixed-tmixed*smixed) c + /(dsqrt(ttmixed-tmixed**2)*dsqrt(ssmixed- c + smixed**2)) c c Calculates new properties in each layer c avg1t=Pconv*tmixed+Pnoconv*til1t avg1s=Pconv*smixed+Pnoconv*til1s sigma1t=sqrt(Pconv*ttmixed+Pnoconv*til1tt-avg1t**2) sigma1s=sqrt(Pconv*ssmixed+Pnoconv*til1ss-avg1s**2) cs1=(Pconv*tsmixed+Pnoconv*til1ts-avg1t*avg1s)/ + (sigma1t*sigma1s) avg2t=Pconv*tmixed+Pnoconv*til2t avg2s=Pconv*smixed+Pnoconv*til2s sigma2t=sqrt(Pconv*ttmixed+Pnoconv*til2tt-avg2t**2) sigma2s=sqrt(Pconv*ssmixed+Pnoconv*til2ss-avg2s**2) cs2=(Pconv*tsmixed+Pnoconv*til2ts-avg2t*avg2s)/ + (sigma2t*sigma2s) avgTMP(il)=avg1t avgTMP(il+1)=avg2t avgS(il)=avg1s avgS(il+1)=avg2s sigmaTMP(il)=sigma1t sigmaS(il)=sigma1s sigmaTMP(il+1)=sigma2t sigmaS(il+1)=sigma2s corr(il)=cs1 corr(il+1)=cs2 goto 10 endif !end mixing enddo !il if (icheck.eq.0) goto 20 enddo 20 continue return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE POTMP(PRESS,TEMP,S,RP,POTEMP) C C TITLE: C ***** C C POTMP -- CALCULATE POTENTIAL TEMPERATURE FOR AN ARBITRARY C REFERENCE PRESSURE C C PURPOSE: C ******* C C TO CALCULATE POTENTIAL TEMPERATURE C C REF: N.P. FOFONOFF C DEEP SEA RESEARCH C IN PRESS NOV 1976 C C PARAMETERS: C ********** C C PRESS -> PRESSURE IN DECIBARS C TEMP -> TEMPERATURE IN CELSIUS DEGREES C S -> SALINITY PSS 78 C RP -> REFERENCE PRESSURE IN DECIBARS C (0.0 FOR THE QUANTITY THETA) C POTEMP <- POTENTIAL TEMPERATURE (DEG C) C REAL*8 PRESS,TEMP,S,RP,POTEMP C C VARIABLES: C ********* C INTEGER I,J,N REAL*4 DP,P,Q,R1,R2,R3,R4,R5,S1,T,X C C CODE: C **** C S1 = S-35.0 P = PRESS T = TEMP C DP = RP - P N = IFIX(ABS(DP)/1000.) + 1 DP = DP/FLOAT(N) C DO 10 I=1,N DO 20 J=1,4 C R1 = ((-2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P R2 = (2.7759E-12*T-1.1351E-10)*S1 R3 = ((-5.4481E-14*T+8.733E-12)*T-6.7795E-10)*T R4 = (R1+(R2+R3+1.8741E-8))*P+(-4.2393E-8*T+1.8932E-6)*S1 R5 = R4+((6.6228E-10*T-6.836E-8)*T+8.5258E-6)*T+3.5803E-5 C X = DP*R5 C GO TO (100,200,300,400),J C 100 CONTINUE T = T+.5*X Q = X P = P + .5*DP GO TO 20 C 200 CONTINUE T = T + .29298322*(X-Q) Q = .58578644*X + .121320344*Q GO TO 20 C 300 CONTINUE T = T + 1.707106781*(X-Q) Q = 3.414213562*X - 4.121320344*Q P = P + .5*DP GO TO 20 C 400 CONTINUE T = T + (X-2.0*Q)/6.0 20 CONTINUE 10 CONTINUE C POTEMP = T RETURN C C END POTMP C END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE RHO(REFPRS,PRSDB,TEMP,SALTY,RHOSTP) C C TITLE: C ***** C C RHO -- CALCULATE DENSITY USING INTERNATIONAL EQUATION OF C STATE C C FROM TEXT FURNISHED BY J. GIESKES C C PARAMETERS: C ********** C C REFPRS -> REFERENCE PRESSURE IN DECIBARS C PRSDB -> INSITU PRESSURE IN DECIBARS C TEMP -> TEMPERATURE IN CELSIUS DEGREES C SALTY -> SALINITY PSS 78 C RHOSTP <- DENSITY IN Kg/M**3 : RHO(S,T,P) C REAL*8 REFPRS,PRSDB,TEMP,SALTY,RHOSTP C C VARIABLES: C ********* C REAL*8 BARS,TERMA,TERMB,KST0,KW,RHOW,KSTP C C CODE: C **** C BARS = REFPRS*.10 C /* REFERENCE PRESSURE IN BARS */ C C /* CALCULATE DENSITY AT GIVEN SALINITY, TEMP AND PRESSURE=0.0 */ C RHOW = -.00909529+TEMP*(1.001685E-4+TEMP*(-1.120083E-6+ 1 TEMP*6.536332E-9)) RHOW = 999.842594+TEMP*(.06793952+TEMP*RHOW) C /* RHOW = DENSITY OF PURE WATER Kg/M**3 */ C KW = TEMP*(-.0040899+TEMP*(7.6438E-5+TEMP*(-8.2467E-7 1 +TEMP*5.3875E-9))) C /* PURE WATER SECANT BULK MODULUS */ C KW = (KW+.824493)*SALTY C KST0 = (-.00572466+TEMP*(1.0227E-4+TEMP*(-1.6546E-6))) 1 * ABS(SALTY)**1.5 C /* K(S,T,0) */ C RHOSTP = RHOW + KW + KST0 + 4.8314E-4*SALTY*SALTY C /* RHO(S,T,P) Kg/M**3 */ C C SELECT /* IF REFERENCE PRESSURE=0.0, DONE */ C (REFERENCE PRESSURE = 0.0): IF(REFPRS.EQ.0.0) GO TO 999 C (OTHERWISE): /* CALCULATE PRESSURE EFFECT */ C C /* RHO(S,T,0)/(1.0-P/K(S,T,P)) */ C KW = TEMP*(148.4206+TEMP*(-2.327105+TEMP 1 *(.01360477+TEMP*(-5.155288E-5)))) KW = KW + 19652.21 C KST0 = (54.6746+TEMP*(-.603459+TEMP*(.0109987+TEMP 1 *(-6.167E-5))))*SALTY C KST0 = KST0 + KW+(.07944+TEMP*(.016483+TEMP*(-5.3009E-4))) 1 *ABS(SALTY)**1.5 C C /* CALCULATE PRESSURE TERMS */ C TERMA = 3.239908+TEMP*(.00143713+TEMP*(1.16092E-4+TEMP 1 *(-5.77905E-7))) TERMA = TERMA + (.0022838+TEMP*(-1.0981E-5+TEMP 1 *(-1.6078E-6)))*SALTY TERMA = TERMA + 1.91075E-4*ABS(SALTY)**1.5 C TERMB = 8.50935E-5+TEMP*(-6.12293E-6+TEMP*5.2787E-8) TERMB = TERMB + (-9.9348E-7+TEMP*(2.0816E-8+TEMP 1 * 9.1697E-10))*SALTY C KSTP = KST0+BARS*(TERMA + BARS*TERMB) C /* SECANT BULK MODULUS K(S,T,P) */ C RHOSTP = RHOSTP/(1.0-BARS/KSTP) C C END SELECT C 999 CONTINUE RETURN C C END RHO C END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE RHO_LIN(t,s,k,rho) integer k,NL real*8 p,t,s,rho,p1 parameter (NL=1000) real*8 alfa(NL),beta(NL),t0(NL),s0(NL) common /leos/alfa,beta,t0,s0 rho=( -(t-t0(k))*alfa(k) + (s-s0(k))*beta(k) ) RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc