! ! AC 2011-01-05 ! SUBROUTINE calltherm_interface (firstcall, & & long,lati,zzlev,zzlay, & & ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, & & pplay,pplev,pphi,zpopsk, & & pdu_th,pdv_th,pdt_th,pdq_th,lmax,zmaxth,pbl_dtke,hfmax,wstar) USE ioipsl_getincom implicit none #include "callkeys.h" #include "dimensions.h" #include "dimphys.h" #include "comcstfi.h" #include "tracer.h" !-------------------------------------------------------- ! Input Variables !-------------------------------------------------------- REAL, INTENT(IN) :: ptimestep REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1),pplay(ngridmx,nlayermx) REAL, INTENT(IN) :: pphi(ngridmx,nlayermx) REAL, INTENT(IN) :: pu(ngridmx,nlayermx),pv(ngridmx,nlayermx) REAL, INTENT(IN) :: pt(ngridmx,nlayermx),pq(ngridmx,nlayermx,nqmx) REAL, INTENT(IN) :: zzlay(ngridmx,nlayermx) REAL, INTENT(IN) :: zzlev(ngridmx,nlayermx+1) LOGICAL, INTENT(IN) :: firstcall REAL, INTENT(IN) :: pdu(ngridmx,nlayermx),pdv(ngridmx,nlayermx) REAL, INTENT(IN) :: pdq(ngridmx,nlayermx,nqmx),pdt(ngridmx,nlayermx) REAL, INTENT(IN) :: q2(ngridmx,nlayermx+1) REAL, INTENT(IN) :: long(ngridmx),lati(ngridmx) REAL, INTENT(IN) :: zpopsk(ngridmx,nlayermx) !-------------------------------------------------------- ! Output Variables !-------------------------------------------------------- REAL, INTENT(OUT) :: pdu_th(ngridmx,nlayermx) REAL, INTENT(OUT) :: pdv_th(ngridmx,nlayermx) REAL, INTENT(OUT) :: pdt_th(ngridmx,nlayermx) REAL, INTENT(OUT) :: pdq_th(ngridmx,nlayermx,nqmx) INTEGER, INTENT(OUT) :: lmax(ngridmx) REAL, INTENT(OUT) :: zmaxth(ngridmx) REAL, INTENT(OUT) :: pbl_dtke(ngridmx,nlayermx+1) REAL, INTENT(OUT) :: wstar(ngridmx) !-------------------------------------------------------- ! Thermals local variables !-------------------------------------------------------- REAL zu(ngridmx,nlayermx), zv(ngridmx,nlayermx) REAL zt(ngridmx,nlayermx) REAL d_t_ajs(ngridmx,nlayermx) REAL d_u_ajs(ngridmx,nlayermx), d_q_ajs(ngridmx,nlayermx,nqmx) REAL d_v_ajs(ngridmx,nlayermx) REAL fm_therm(ngridmx,nlayermx+1), entr_therm(ngridmx,nlayermx) REAL detr_therm(ngridmx,nlayermx) REAL zw2(ngridmx,nlayermx+1) REAL fraca(ngridmx,nlayermx+1),zfraca(ngridmx,nlayermx+1) REAL ztla(ngridmx,nlayermx) REAL q_therm(ngridmx,nlayermx), pq_therm(ngridmx,nlayermx,nqmx) REAL q2_therm(ngridmx,nlayermx), dq2_therm(ngridmx,nlayermx) REAL lmax_real(ngridmx) REAL masse(ngridmx,nlayermx) REAL zdz(ngridmx,nlayermx) LOGICAL qtransport_thermals,dtke_thermals INTEGER l,ig,iq,ii(1) CHARACTER (LEN=20) :: modname !-------------------------------------------------------- ! Local variables for sub-timestep !-------------------------------------------------------- REAL d_t_the(ngridmx,nlayermx), d_q_the(ngridmx,nlayermx,nqmx) REAL d_u_the(ngridmx,nlayermx),d_v_the(ngridmx,nlayermx) REAL dq2_the(ngridmx,nlayermx) INTEGER isplit,nsplit_thermals REAL r_aspect_thermals REAL fact REAL zfm_therm(ngridmx,nlayermx+1),zdt REAL zentr_therm(ngridmx,nlayermx),zdetr_therm(ngridmx,nlayermx) REAL zheatFlux(ngridmx,nlayermx) REAL zheatFlux_down(ngridmx,nlayermx) REAL zbuoyancyOut(ngridmx,nlayermx) REAL zbuoyancyEst(ngridmx,nlayermx) REAL zzw2(ngridmx,nlayermx+1) REAL zmax(ngridmx) !-------------------------------------------------------- ! Diagnostics !-------------------------------------------------------- REAL heatFlux(ngridmx,nlayermx) REAL heatFlux_down(ngridmx,nlayermx) REAL buoyancyOut(ngridmx,nlayermx) REAL buoyancyEst(ngridmx,nlayermx) REAL hfmax(ngridmx),wmax(ngridmx) REAL pbl_teta(ngridmx),dteta(ngridmx,nlayermx) !-------------------------------------------------------- ! Theta_m !-------------------------------------------------------- INTEGER ico2 SAVE ico2 ! ********************************************************************** ! Initialization ! ********************************************************************** lmax(:)=0. pdu_th(:,:)=0. pdv_th(:,:)=0. pdt_th(:,:)=0. entr_therm(:,:)=0. detr_therm(:,:)=0. q2_therm(:,:)=0. dq2_therm(:,:)=0. ztla(:,:)=0. pbl_dtke(:,:)=0. fm_therm(:,:)=0. zw2(:,:)=0. fraca(:,:)=0. zfraca(:,:)=0. if (tracer) then pdq_th(:,:,:)=0. end if d_t_ajs(:,:)=0. d_u_ajs(:,:)=0. d_v_ajs(:,:)=0. d_q_ajs(:,:,:)=0. heatFlux(:,:)=0. heatFlux_down(:,:)=0. buoyancyOut(:,:)=0. buoyancyEst(:,:)=0. zmaxth(:)=0. lmax_real(:)=0. ! ********************************************************************** ! Preparing inputs for the thermals ! ********************************************************************** zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep pq_therm(:,:,:)=0. qtransport_thermals=.true. !! default setting !call getin("qtransport_thermals",qtransport_thermals) if(qtransport_thermals) then if(tracer) then pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep endif endif ! dtke_thermals=.false. !! default setting ! !call getin("dtke_thermals",dtke_thermals) ! ! IF(dtke_thermals) THEN ! DO l=1,nlayermx ! q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1)) ! ENDDO ! ENDIF ! ********************************************************************** ! Polar night mixing : theta_m ! ********************************************************************** if(firstcall) then ico2=0 if (tracer) then ! Prepare Special treatment if one of the tracers is CO2 gas do iq=1,nqmx if (noms(iq).eq."co2") then ico2=iq end if enddo endif endif !of if firstcall ! ********************************************************************** ! ********************************************************************** ! ********************************************************************** ! CALLTHERM ! ********************************************************************** ! ********************************************************************** ! ********************************************************************** ! r_aspect_thermals ! ultimately conrols the amount of mass going through the thermals ! decreasing it increases the thermals effect. Tests at gcm resolution ! shows that too low values destabilize the model ! when changing this value, one should check that the surface layer model ! outputs the correct Cd*u and Ch*u through changing the gustiness coefficient beta #ifdef MESOSCALE !! valid for timesteps < 200s nsplit_thermals=2 r_aspect_thermals=0.7 #else nsplit_thermals=35 r_aspect_thermals=1.5 #endif call getin("nsplit_thermals",nsplit_thermals) call getin("r_aspect_thermals",r_aspect_thermals) ! ********************************************************************** ! SUB-TIMESTEP LOOP ! ********************************************************************** zdt=ptimestep/REAL(nsplit_thermals) DO isplit=1,nsplit_thermals ! Initialization of intermediary variables zfm_therm(:,:)=0. zentr_therm(:,:)=0. zdetr_therm(:,:)=0. zheatFlux(:,:)=0. zheatFlux_down(:,:)=0. zbuoyancyOut(:,:)=0. zbuoyancyEst(:,:)=0. zzw2(:,:)=0. zmax(:)=0. lmax(:)=0. d_t_the(:,:)=0. d_u_the(:,:)=0. d_v_the(:,:)=0. dq2_the(:,:)=0. if (nqmx .ne. 0) then d_q_the(:,:,:)=0. endif CALL thermcell_main_mars(zdt & & ,pplay,pplev,pphi,zzlev,zzlay & & ,zu,zv,zt,pq_therm,q2_therm & & ,d_u_the,d_v_the,d_t_the,d_q_the,dq2_the & & ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax & & ,r_aspect_thermals & & ,zzw2,fraca,zpopsk & & ,ztla,zheatFlux,zheatFlux_down & & ,zbuoyancyOut,zbuoyancyEst) fact=1./REAL(nsplit_thermals) d_t_the(:,:)=d_t_the(:,:)*ptimestep*fact ! d_u_the(:,:)=d_u_the(:,:)*fact ! d_v_the(:,:)=d_v_the(:,:)*fact ! dq2_the(:,:)=dq2_the(:,:)*fact if (ico2 .ne. 0) then d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*fact endif zmaxth(:)=zmaxth(:)+zmax(:)*fact lmax_real(:)=lmax_real(:)+float(lmax(:))*fact fm_therm(:,:)=fm_therm(:,:) & & +zfm_therm(:,:)*fact entr_therm(:,:)=entr_therm(:,:) & & +zentr_therm(:,:)*fact detr_therm(:,:)=detr_therm(:,:) & & +zdetr_therm(:,:)*fact zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact heatFlux(:,:)=heatFlux(:,:) & & +zheatFlux(:,:)*fact heatFlux_down(:,:)=heatFlux_down(:,:) & & +zheatFlux_down(:,:)*fact buoyancyOut(:,:)=buoyancyOut(:,:) & & +zbuoyancyOut(:,:)*fact buoyancyEst(:,:)=buoyancyEst(:,:) & & +zbuoyancyEst(:,:)*fact zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact ! accumulation de la tendance d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:) ! d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:) ! d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:) if (ico2 .ne. 0) then d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2) endif ! dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:) ! incrementation des variables meteo zt(:,:) = zt(:,:) + d_t_the(:,:) ! zu(:,:) = zu(:,:) + d_u_the(:,:) ! zv(:,:) = zv(:,:) + d_v_the(:,:) if (ico2 .ne. 0) then pq_therm(:,:,ico2) = & & pq_therm(:,:,ico2) + d_q_the(:,:,ico2)*ptimestep endif ! q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:) ENDDO ! isplit !**************************************************************** ! Now that we have computed total entrainment and detrainment, we can ! advect u, v, and q in thermals. (theta already advected). We can do ! that separatly because u,v,and q are not used in thermcell_main for ! any thermals-related computation : they are purely passive. ! mass of cells do l=1,nlayermx masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g enddo ! thickness of layers do l=1,nlayermx zdz(:,l)=zzlev(:,l+1)-zzlev(:,l) enddo modname='momentum' call thermcell_dqup(ngridmx,nlayermx,ptimestep & & ,fm_therm,entr_therm,detr_therm, & & masse,zu,d_u_ajs,modname,zdz) call thermcell_dqup(ngridmx,nlayermx,ptimestep & & ,fm_therm,entr_therm,detr_therm, & & masse,zv,d_v_ajs,modname,zdz) if (nqmx .ne. 0.) then modname='tracer' DO iq=1,nqmx if (iq .ne. ico2) then call thermcell_dqup(ngridmx,nlayermx,ptimestep & & ,fm_therm,entr_therm,detr_therm, & & masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),modname,zdz) endif ENDDO endif DO ig=1,ngridmx hfmax(ig)=MAXVAL(heatFlux(ig,:)+heatFlux_down(ig,:)) wmax(ig)=MAXVAL(zw2(ig,:)) ENDDO lmax(:)=nint(lmax_real(:)) ! ********************************************************************** ! ********************************************************************** ! ********************************************************************** ! CALLTHERM END ! ********************************************************************** ! ********************************************************************** ! ********************************************************************** ! ********************************************************************** ! Preparing outputs ! ********************************************************************** ! Winds and tracers PDU, PDV, and PDQ are DERIVATIVES pdu_th(:,:)=d_u_ajs(:,:) pdv_th(:,:)=d_v_ajs(:,:) if(qtransport_thermals) then if(tracer) then pdq_th(:,:,:)=d_q_ajs(:,:,:) endif endif ! IF(dtke_thermals) THEN ! DO l=2,nlayermx ! pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l)) ! ENDDO ! ! pbl_dtke(:,1)=0.5*dq2_therm(:,1) ! pbl_dtke(:,nlayermx+1)=0. ! ENDIF ! Temperature PDT is a TENDANCY pdt_th(:,:)=d_t_ajs(:,:)/ptimestep ! ********************************************************************** ! Compute the free convection velocity scale for vdifc ! ********************************************************************** ! Potential temperature gradient dteta(:,nlayermx)=0. DO l=1,nlayermx-1 DO ig=1, ngridmx dteta(ig,l) = ((zt(ig,l+1)-zt(ig,l))/zpopsk(ig,l)) & & /(zzlay(ig,l+1)-zzlay(ig,l)) ENDDO ENDDO ! Computation of the pbl mixed layer temperature DO ig=1, ngridmx ii=MINLOC(abs(dteta(ig,1:lmax(ig)))) pbl_teta(ig) = zt(ig,ii(1))/zpopsk(ig,ii(1)) ENDDO ! We follow Spiga et. al 2010 (QJRMS) ! ------------ DO ig=1, ngridmx IF (zmax(ig) .gt. 0.) THEN wstar(ig)=(g*zmaxth(ig)*hfmax(ig)/pbl_teta(ig))**(1./3.) ELSE wstar(ig)=0. ENDIF ENDDO ! ********************************************************************** ! Diagnostics ! ********************************************************************** if(outptherm) then if (ngridmx .eq. 1) then call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',& & 'kg/m-2',1,entr_therm) call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',& & 'kg/m-2',1,detr_therm) call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',& & 'kg/m-2',1,fm_therm) call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',& & 'm/s',1,zw2) call WRITEDIAGFI(ngridmx,'heatFlux_up','heatFlux_updraft',& & 'SI',1,heatFlux) call WRITEDIAGFI(ngridmx,'heatFlux_down','heatFlux_downdraft',& & 'SI',1,heatFlux_down) call WRITEDIAGFI(ngridmx,'fraca','fraction coverage',& & 'percent',1,fraca) call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',& & 'm.s-2',1,buoyancyOut) call WRITEDIAGFI(ngridmx,'buoyancyEst','buoyancyEst',& & 'm.s-2',1,buoyancyEst) call WRITEDIAGFI(ngridmx,'d_t_th', & & 'tendance temp TH','K',1,d_t_ajs) call WRITEDIAGFI(ngridmx,'zmax', & & 'pbl height','m',0,zmaxth) else call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',& & 'kg/m-2',3,entr_therm) call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',& & 'kg/m-2',3,detr_therm) call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',& & 'kg/m-2',3,fm_therm) call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',& & 'm/s',3,zw2) call WRITEDIAGFI(ngridmx,'heatFlux','heatFlux',& & 'SI',3,heatFlux) call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',& & 'SI',3,buoyancyOut) call WRITEDIAGFI(ngridmx,'d_t_th', & & 'tendance temp TH','K',3,d_t_ajs) endif endif END