subroutine moistadj(t, pq, pplev, pplay, dtmana, dqmana, ptimestep, rneb) use watercommon_h, only: To, RLVTT, RCPD implicit none !===================================================================== ! ! Purpose ! ------- ! Calculates moist convective adjustment by the method of Manabe. ! ! Authors ! ------- ! Adapted from the LMDTERRE code by R. Wordsworth (2010) ! Original author Z. X. Li (1993) ! !===================================================================== #include "dimensions.h" #include "dimphys.h" #include "tracer.h" #include "comcstfi.h" ! Pre-arguments (for universal model) real pq(ngridmx,nlayermx,nqmx) ! tracer (kg/kg) REAL pdq(ngridmx,nlayermx,nqmx) real dqmana(ngridmx,nlayermx,nqmx) ! tendency of tracers (kg/kg.s-1) REAL dtmana(ngridmx,nlayermx) ! temperature increment ! Arguments REAL t(ngridmx,nlayermx) ! temperature (K) REAL q(ngridmx,nlayermx) ! humidite specifique (kg/kg) REAL pplev(ngridmx,nlayermx+1) ! pression a inter-couche (Pa) REAL pplay(ngridmx,nlayermx) ! pression au milieu de couche (Pa) REAL d_t(ngridmx,nlayermx) ! temperature increment REAL d_q(ngridmx,nlayermx) ! incrementation pour vapeur d'eau REAL d_ql(ngridmx,nlayermx) ! incrementation pour l'eau liquide REAL rneb(ngridmx,nlayermx) ! cloud fraction REAL ptimestep ! REAL t_coup ! PARAMETER (t_coup=234.0) REAL seuil_vap PARAMETER (seuil_vap=1.0E-10) ! Local variables INTEGER i, k, iq INTEGER k1, k1p, k2, k2p LOGICAL itest(ngridmx) REAL delta_q(ngridmx, nlayermx) REAL cp_new_t(nlayermx) REAL cp_delta_t(nlayermx) REAL new_qb(nlayermx) REAL v_cptj(nlayermx), v_cptjk1, v_ssig REAL v_cptt(ngridmx,nlayermx), v_p, v_t REAL v_qs(ngridmx,nlayermx), v_qsd(ngridmx,nlayermx) REAL zq1(ngridmx), zq2(ngridmx) REAL gamcpdz(ngridmx,2:nlayermx) REAL zdp, zdpm REAL zsat ! super-saturation REAL zflo ! flotabilite REAL local_q(ngridmx,nlayermx),local_t(ngridmx,nlayermx) REAL zdelta, zcor, zcvm5 REAL dEtot, dqtot, masse ! conservation diagnostics real dL1tot, dL2tot ! Indices of water vapour and water ice tracers INTEGER,SAVE :: i_h2o=0 ! water vapour INTEGER,SAVE :: i_ice=0 ! water ice LOGICAL firstcall SAVE firstcall DATA firstcall /.TRUE./ IF (firstcall) THEN i_h2o=igcm_h2o_vap i_ice=igcm_h2o_ice write(*,*) "rain: i_ice=",i_ice write(*,*) " i_h2o=",i_h2o firstcall = .FALSE. ENDIF ! GCM -----> subroutine variables DO k = 1, nlayermx DO i = 1, ngridmx q(i,k) = pq(i,k,i_h2o) if(q(i,k).lt.0.)then q(i,k)=0.0 endif DO iq = 1, nqmx dqmana(i,k,iq)=0.0 ENDDO ENDDO ENDDO DO k = 1, nlayermx DO i = 1, ngridmx local_q(i,k) = q(i,k) local_t(i,k) = t(i,k) rneb(i,k) = 0.0 d_ql(i,k) = 0.0 d_t(i,k) = 0.0 d_q(i,k) = 0.0 ENDDO new_qb(k)=0.0 ENDDO ! Calculate v_cptt DO k = 1, nlayermx DO i = 1, ngridmx v_cptt(i,k) = RCPD * local_t(i,k) ENDDO ENDDO DO k = 1, nlayermx DO i = 1, ngridmx v_cptt(i,k) = RCPD * local_t(i,k) v_t = local_t(i,k) v_p = pplay(i,k) call watersat(v_t,v_p,v_qs(i,k)) call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k)) ENDDO ENDDO ! TEST: RH DIAGNOSTIC ! DO k = 1, nlayermx ! DO i = 1, ngridmx ! v_t = local_t(i,k) ! IF (v_t.LT.To) THEN ! print*,'RHs=',q(i,k) / v_qs(i,k) ! ELSE ! print*,'RHl=',q(i,k) / v_qs(i,k) ! ENDIF ! ENDDO ! ENDDO ! Calculate Gamma * Cp * dz: (gamma is the critical gradient) DO k = 2, nlayermx DO i = 1, ngridmx zdp = pplev(i,k)-pplev(i,k+1) zdpm = pplev(i,k-1)-pplev(i,k) ! gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) * gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * & (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp) & +RLVTT /(zdpm+zdp) * & (v_qs(i,k-1)*zdpm + v_qs(i,k)*zdp) & )* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) & / (1.0+(v_qsd(i,k-1)*zdpm+ & v_qsd(i,k)*zdp)/(zdpm+zdp) ) ENDDO ENDDO !------------------------------------ modification of unstable profile DO 9999 i = 1, ngridmx itest(i) = .FALSE. ! print*,'we in the loop' ! stop k1 = 0 k2 = 1 810 CONTINUE ! look for k1, the base of the column k2 = k2 + 1 IF (k2 .GT. nlayermx) GOTO 9999 zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2) zsat=(local_q(i,k2-1)-v_qs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2)) & +(local_q(i,k2)-v_qs(i,k2))*(pplev(i,k2)-pplev(i,k2+1)) IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 ) GOTO 810 k1 = k2 - 1 itest(i) = .TRUE. 820 CONTINUE !! look for k2, the top of the column IF (k2 .EQ. nlayermx) GOTO 821 k2p = k2 + 1 zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-v_qs(i,k2p)) zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p) IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 821 k2 = k2p GOTO 820 821 CONTINUE !------------------------------------------------------ local adjustment 830 CONTINUE ! actual adjustment v_cptj(k1) = 0.0 zdp = pplev(i,k1)-pplev(i,k1+1) v_cptjk1 = ( (1.0+v_qsd(i,k1))*(v_cptt(i,k1)+v_cptj(k1)) & + RLVTT*(local_q(i,k1)-v_qs(i,k1)) ) * zdp v_ssig = zdp * (1.0+v_qsd(i,k1)) k1p = k1 + 1 DO k = k1p, k2 zdp = pplev(i,k)-pplev(i,k+1) v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k) v_cptjk1 = v_cptjk1 + zdp & * ( (1.0+v_qsd(i, k))*(v_cptt(i,k)+v_cptj(k)) & + RLVTT*(local_q(i,k)-v_qs(i,k)) ) v_ssig = v_ssig + zdp *(1.0+v_qsd(i,k)) ENDDO ! this right here is where the adjustment is done??? DO k = k1, k2 cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k) cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k) new_qb(k) = v_qs(i,k) + v_qsd(i,k)*cp_delta_t(k)/RLVTT local_q(i,k) = new_qb(k) local_t(i,k) = cp_new_t(k) / RCPD ! print*,'v_qs in loop=',v_qs ! print*,'v_qsd in loop=',v_qsd ! print*,'new_qb in loop=',new_qb ! print*,'cp_delta_t in loop=',cp_delta_t ENDDO !--------------------------------------------------- sounding downwards ! -- we refine the prognostic variables in ! -- the layer about to be adjusted DO k = k1, k2 v_cptt(i,k) = RCPD * local_t(i,k) v_t = local_t(i,k) v_p = pplay(i,k) ! IF (v_t.LT.t_coup) THEN ! v_qs(i,k) = qsats(v_t) / v_p ! v_qsd(i,k) = dqsats(v_t,v_qs(i,k)) ! ELSE ! v_qs(i,k) = qsatl(v_t) / v_p ! v_qsd(i,k) = dqsatl(v_t,v_qs(i,k)) ! ENDIF call watersat(v_t,v_p,v_qs(i,k)) call watersat_grad(v_t,v_qs(i,k),v_qsd(i,k)) ENDDO DO k = 2, nlayermx zdpm = pplev(i,k-1) - pplev(i,k) zdp = pplev(i,k) - pplev(i,k+1) ! gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) * gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * & (v_cptt(i,k-1)*zdpm+v_cptt(i,k)*zdp) & +RLVTT /(zdpm+zdp) * & (v_qs(i,k-1)*zdpm+v_qs(i,k)*zdp) & )* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) & / (1.0+(v_qsd(i,k-1)*zdpm+v_qsd(i,k)*zdp) & /(zdpm+zdp) ) ENDDO ! Test to see if we've reached the bottom IF (k1 .EQ. 1) GOTO 841 ! yes we have! zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1) zsat=(local_q(i,k1-1)-v_qs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1)) & + (local_q(i,k1)-v_qs(i,k1))*(pplev(i,k1)-pplev(i,k1+1)) IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have! 840 CONTINUE k1 = k1 - 1 IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995) zsat = zsat + (local_q(i,k1-1)-v_qs(i,k1-1)) & *(pplev(i,k1-1)-pplev(i,k1)) zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1) IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN GOTO 840 ELSE GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995) ENDIF 841 CONTINUE GOTO 810 ! look for other layers higher up 9999 CONTINUE ! loop over all the points ! print*,'k1=',k1 ! print*,'k2=',k2 ! print*,'local_t=',local_t ! print*,'v_cptt=',v_cptt ! print*,'gamcpdz=',gamcpdz !----------------------------------------------------------------------- ! Determine the cloud fraction (hypothese: la nebulosite a lieu ! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement): DO k = 1, nlayermx DO i = 1, ngridmx IF (itest(i)) THEN delta_q(i,k) = local_q(i,k) - q(i,k) IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0 ENDIF ENDDO ENDDO ! Distribuer l'eau condensee en eau liquide nuageuse (hypothese: ! l'eau liquide est distribuee aux endroits ou la vapeur d'eau ! diminue et d'une maniere proportionnelle a cet diminution): DO i = 1, ngridmx IF (itest(i)) THEN zq1(i) = 0.0 zq2(i) = 0.0 ENDIF ENDDO DO k = 1, nlayermx DO i = 1, ngridmx IF (itest(i)) THEN zdp = pplev(i,k)-pplev(i,k+1) zq1(i) = zq1(i) - delta_q(i,k) * zdp zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp ENDIF ENDDO ENDDO DO k = 1, nlayermx DO i = 1, ngridmx IF (itest(i)) THEN IF (zq2(i).NE.0.0) & d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i) ENDIF ENDDO ENDDO ! print*,'local_q BEFORE=',local_q DO k = 1, nlayermx DO i = 1, ngridmx local_q(i, k) = MAX(local_q(i, k), seuil_vap) ENDDO ENDDO DO k = 1, nlayermx DO i = 1, ngridmx d_t(i,k) = local_t(i,k) - t(i,k) d_q(i,k) = local_q(i,k) - q(i,k) ENDDO ENDDO ! now subroutine -----> GCM variables DO k = 1, nlayermx DO i = 1, ngridmx dtmana(i,k) = d_t(i,k)/ptimestep dqmana(i,k,i_h2o) = d_q(i,k)/ptimestep dqmana(i,k,i_ice) = d_ql(i,k)/ptimestep ENDDO ENDDO ! print*,'IN MANABE:' ! print*,'pplev=',pplev ! print*,'t=',t ! print*,'d_t=',d_t ! print*,'d_q=',d_q ! print*,'local_q=',local_q ! print*,'q=',q ! print*,'v_qs(i,k)=',v_qs ! print*,'v_qsd(i,k)=',v_qsd ! print*,'cp_delta_t(k)=',cp_delta_t ! print*,'d_ql=',d_ql ! print*,'delta_q=',delta_q ! print*,'zq1=',zq1 ! print*,'zq2=',zq2 !! print*,'i_h2o=',i_h2o ! print*,'i_ice=',i_ice ! ! print*,'IN MANABE:' ! print*,'d_q=',d_q ! print*,'d_ql=',d_ql ! print*,'dtmana=',d_t ! stop ! print*,'gamcpdz at end=',gamcpdz ! stop ! Some conservation diagnostics... ! dEtot=0.0 ! dL1tot=0.0 ! dL2tot=0.0 ! dqtot=0.0 ! masse=0.0 ! DO k = 1, nlayermx ! DO i = 1, ngridmx ! ! masse = (pplev(i,k) - pplev(i,k+1))/g ! ! dEtot = dEtot + cpp*d_t(i,k)*masse ! dL1tot = dL1tot + RLVTT*d_ql(i,k)*masse ! dL2tot = dL2tot + RLVTT*d_q(i,k)*masse ! is this line necessary? ! ! dqtot = dqtot + (d_q(i,k) + d_ql(i,k))*masse ! ! ENDDO ! ENDDO ! print*,'In manabe energy change=',dEtot ! print*,'In manabe condense energy change 1 =',dL1tot ! print*,'In manabe condense energy change 2 =',dL2tot ! print*,'In manabe water change=',dqtot RETURN END