! ! $Header$ ! SUBROUTINE vlspltgen_loc( q,pente_max,masse,w,pbaru,pbarv, & pdt, p,pk,teta ) ! ! Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron ! ! ******************************************************************** ! Schema d'advection " pseudo amont " . ! + test sur humidite specifique: Q advecte< Qsat aval ! (F. Codron, 10/99) ! ******************************************************************** ! q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... ! ! pente_max facteur de limitation des pentes: 2 en general ! 0 pour un schema amont ! pbaru,pbarv,w flux de masse en u ,v ,w ! pdt pas de temps ! ! teta temperature potentielle, p pression aux interfaces, ! pk exner au milieu des couches necessaire pour calculer Qsat ! -------------------------------------------------------------------- USE parallel_lmdz USE mod_hallo USE Write_Field_loc USE VAMPIR ! ! CRisi: on rajoute variables utiles d'infotrac USE infotrac, ONLY : nqtot, tracers, isoCheck USE vlspltgen_mod USE comconst_mod, ONLY: cpp USE logic_mod, ONLY: adv_qsat_liq IMPLICIT NONE ! include "dimensions.h" include "paramet.h" ! ! Arguments: ! ---------- REAL :: masse(ijb_u:ije_u,llm),pente_max REAL :: pbaru( ijb_u:ije_u,llm ),pbarv( ijb_v:ije_v,llm) REAL :: q(ijb_u:ije_u,llm,nqtot) REAL :: w(ijb_u:ije_u,llm),pdt REAL :: p(ijb_u:ije_u,llmp1),teta(ijb_u:ije_u,llm) REAL :: pk(ijb_u:ije_u,llm) ! ! Local ! --------- ! INTEGER :: ij,l ! REAL :: zzpbar, zzw REAL :: qmin,qmax DATA qmin,qmax/0.,1.e33/ !--pour rapport de melange saturant-- REAL :: rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play REAL :: ptarg,pdelarg,foeew,zdelta REAL :: tempe(ijb_u:ije_u) INTEGER :: ijb,ije,iq,iq2,ifils LOGICAL, SAVE :: firstcall=.TRUE. !$OMP THREADPRIVATE(firstcall) type(request),SAVE :: MyRequest1 !$OMP THREADPRIVATE(MyRequest1) type(request),SAVE :: MyRequest2 !$OMP THREADPRIVATE(MyRequest2) ! fonction psat(T) FOEEW ( PTARG,PDELARG ) = EXP ( & (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) ) r2es = 380.11733 r3les = 17.269 r3ies = 21.875 r4les = 35.86 r4ies = 7.66 retv = 0.6077667 rtt = 273.16 ! Allocate variables depending on dynamic variable nqtot IF (firstcall) THEN firstcall=.FALSE. END IF !-- Calcul de Qsat en chaque point !-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2 ! pour eviter une exponentielle. call SetTag(MyRequest1,100) call SetTag(MyRequest2,101) ijb=ij_begin-iip1 ije=ij_end+iip1 if (pole_nord) ijb=ij_begin if (pole_sud) ije=ij_end !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, llm DO ij = ijb, ije tempe(ij) = teta(ij,l) * pk(ij,l) /cpp ENDDO DO ij = ijb, ije IF (adv_qsat_liq) THEN zdelta = 0. ELSE zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) ) ENDIF play = 0.5*(p(ij,l)+p(ij,l+1)) qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play ) qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) ) ENDDO ENDDO !$OMP END DO NOWAIT ! PRINT*,'Debut vlsplt version debug sans vlyqs' zzpbar = 0.5 * pdt zzw = pdt ijb=ij_begin ije=ij_end if (pole_nord) ijb=ijb+iip1 if (pole_sud) ije=ije-iip1 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm DO ij = ijb,ije mu(ij,l)=pbaru(ij,l) * zzpbar ENDDO ENDDO !$OMP END DO NOWAIT ijb=ij_begin-iip1 ije=ij_end if (pole_nord) ijb=ij_begin if (pole_sud) ije=ij_end-iip1 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm DO ij=ijb,ije mv(ij,l)=pbarv(ij,l) * zzpbar ENDDO ENDDO !$OMP END DO NOWAIT ijb=ij_begin ije=ij_end DO iq=1,nqtot !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm DO ij=ijb,ije mw(ij,l,iq)=w(ij,l) * zzw ENDDO ENDDO !$OMP END DO NOWAIT ENDDO DO iq=1,nqtot !$OMP MASTER DO ij=ijb,ije mw(ij,llm+1,iq)=0. ENDDO !$OMP END MASTER ENDDO ! CALL SCOPY(ijp1llm,q,1,zq,1) ! CALL SCOPY(ijp1llm,masse,1,zm,1) ijb=ij_begin ije=ij_end DO iq=1,nqtot !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm zq(ijb:ije,l,iq)=q(ijb:ije,l,iq) zm(ijb:ije,l,iq)=masse(ijb:ije,l) ENDDO !$OMP END DO NOWAIT ENDDO #ifdef DEBUG_IO CALL WriteField_u('mu',mu) CALL WriteField_v('mv',mv) CALL WriteField_u('mw',mw) CALL WriteField_u('qsat',qsat) #endif ! ! verif temporaire ijb=ij_begin ije=ij_end call check_isotopes(zq,ijb,ije,'vlspltgen_loc 191') !$OMP BARRIER DO iq=1,nqtot ! ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air IF(tracers(iq)%parent /= 'air') CYCLE ! !write(*,*) 'vlspltgen 192: iq,iadv=',iq,tracers(iq)%iadv #ifdef DEBUG_IO CALL WriteField_u('zq',zq(:,:,iq)) CALL WriteField_u('zm',zm(:,:,iq)) #endif SELECT CASE(tracers(iq)%iadv) CASE(0); CYCLE CASE(10) #ifdef _ADV_HALO ! CRisi: on ajoute les nombres de fils et tableaux des fils ! On suppose qu'on ne peut advecter les fils que par le schéma 10. call vlx_loc(zq,pente_max,zm,mu, & ij_begin,ij_begin+2*iip1-1,iq) call vlx_loc(zq,pente_max,zm,mu, & ij_end-2*iip1+1,ij_end,iq) #else call vlx_loc(zq,pente_max,zm,mu, & ij_begin,ij_end,iq) #endif !$OMP MASTER call VTb(VTHallo) !$OMP END MASTER call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1) call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1) ! CRisi do ifils=1,tracers(iq)%nqDescen iq2=tracers(iq)%iqDescen(ifils) call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1) call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1) enddo !$OMP MASTER call VTe(VTHallo) !$OMP END MASTER CASE(14) #ifdef _ADV_HALO call vlxqs_loc(zq,pente_max,zm,mu, & qsat,ij_begin,ij_begin+2*iip1-1,iq) call vlxqs_loc(zq,pente_max,zm,mu, & qsat,ij_end-2*iip1+1,ij_end,iq) #else call vlxqs_loc(zq,pente_max,zm,mu, & qsat,ij_begin,ij_end,iq) #endif !$OMP MASTER call VTb(VTHallo) !$OMP END MASTER call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1) call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1) do ifils=1,tracers(iq)%nqDescen iq2=tracers(iq)%iqDescen(ifils) call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1) call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1) enddo !$OMP MASTER call VTe(VTHallo) !$OMP END MASTER CASE DEFAULT CALL abort_gcm("vlspltgen_loc","schema non parallelise",1) END SELECT enddo !DO iq=1,nqtot !$OMP BARRIER !$OMP MASTER call VTb(VTHallo) !$OMP END MASTER call SendRequest(MyRequest1) !$OMP MASTER call VTe(VTHallo) !$OMP END MASTER !$OMP BARRIER ! ! verif temporaire ijb=ij_begin-2*iip1 ije=ij_end+2*iip1 if (pole_nord) ijb=ij_begin if (pole_sud) ije=ij_end call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 280') do iq=1,nqtot IF(tracers(iq)%parent /= 'air') CYCLE ! !write(*,*) 'vlspltgen 279: iq=',iq SELECT CASE(tracers(iq)%iadv) CASE(0); CYCLE CASE(10) #ifdef _ADV_HALLO call vlx_loc(zq,pente_max,zm,mu, & ij_begin+2*iip1,ij_end-2*iip1,iq) #endif CASE(14) #ifdef _ADV_HALLO call vlxqs_loc(zq,pente_max,zm,mu, & qsat,ij_begin+2*iip1,ij_end-2*iip1,iq) #endif CASE DEFAULT CALL abort_gcm("vlspltgen_p","schema non parallelise",1) END SELECT enddo !$OMP BARRIER !$OMP MASTER call VTb(VTHallo) !$OMP END MASTER ! call WaitRecvRequest(MyRequest1) ! call WaitSendRequest(MyRequest1) !$OMP BARRIER call WaitRequest(MyRequest1) !$OMP MASTER call VTe(VTHallo) !$OMP END MASTER !$OMP BARRIER IF(isoCheck) THEN call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 326') ijb=ij_begin-2*iip1 ije=ij_end+2*iip1 if (pole_nord) ijb=ij_begin if (pole_sud) ije=ij_end call check_isotopes(zq,ijb,ije,'vlspltgen_loc 336') END IF do iq = 1, nqtot IF(tracers(iq)%parent /= 'air') CYCLE ! !write(*,*) 'vlspltgen 321: iq=',iq #ifdef DEBUG_IO CALL WriteField_u('zq',zq(:,:,iq)) CALL WriteField_u('zm',zm(:,:,iq)) #endif SELECT CASE(tracers(iq)%iadv) CASE(0); CYCLE CASE(10); call vly_loc(zq,pente_max,zm,mv, iq) CASE(14); call vlyqs_loc(zq,pente_max,zm,mv,qsat,iq) CASE DEFAULT CALL abort_gcm("vlspltgen_p","schema non parallelise",1) END SELECT enddo call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 357') do iq = 1, nqtot IF(tracers(iq)%parent /= 'air') CYCLE ! !write(*,*) 'vlspltgen 349: iq=',iq #ifdef DEBUG_IO CALL WriteField_u('zq',zq(:,:,iq)) CALL WriteField_u('zm',zm(:,:,iq)) #endif SELECT CASE(tracers(iq)%iadv) CASE(0); CYCLE CASE(10,14) !$OMP BARRIER #ifdef _ADV_HALLO call vlz_loc(zq,pente_max,zm,mw, & ij_begin,ij_begin+2*iip1-1,iq) call vlz_loc(zq,pente_max,zm,mw, & ij_end-2*iip1+1,ij_end,iq) #else call vlz_loc(zq,pente_max,zm,mw, & ij_begin,ij_end,iq) #endif !$OMP BARRIER !$OMP MASTER call VTb(VTHallo) !$OMP END MASTER call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest2) call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest2) ! ! CRisi do ifils=1,tracers(iq)%nqDescen iq2=tracers(iq)%iqDescen(ifils) call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest2) call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest2) enddo !$OMP MASTER call VTe(VTHallo) !$OMP END MASTER !$OMP BARRIER CASE DEFAULT CALL abort_gcm("vlspltgen_p","schema non parallelise",1) END SELECT enddo !$OMP BARRIER !$OMP MASTER call VTb(VTHallo) !$OMP END MASTER call SendRequest(MyRequest2) !$OMP MASTER call VTe(VTHallo) !$OMP END MASTER call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 429') !$OMP BARRIER do iq=1,nqtot IF(tracers(iq)%parent /= 'air') CYCLE ! !write(*,*) 'vlspltgen 409: iq=',iq SELECT CASE(tracers(iq)%iadv) CASE(0); CYCLE CASE(10,14) !$OMP BARRIER #ifdef _ADV_HALLO call vlz_loc(zq,pente_max,zm,mw, & ij_begin+2*iip1,ij_end-2*iip1,iq) #endif !$OMP BARRIER CASE DEFAULT CALL abort_gcm("vlspltgen_p","schema non parallelise",1) END SELECT enddo ! !write(*,*) 'vlspltgen_loc 476' !$OMP BARRIER ! !write(*,*) 'vlspltgen_loc 477' !$OMP MASTER call VTb(VTHallo) !$OMP END MASTER ! call WaitRecvRequest(MyRequest2) ! call WaitSendRequest(MyRequest2) !$OMP BARRIER CALL WaitRequest(MyRequest2) !$OMP MASTER call VTe(VTHallo) !$OMP END MASTER !$OMP BARRIER ! !write(*,*) 'vlspltgen_loc 494' call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 461') do iq=1,nqtot IF(tracers(iq)%parent /= 'air') CYCLE ! !write(*,*) 'vlspltgen 449: iq=',iq #ifdef DEBUG_IO CALL WriteField_u('zq',zq(:,:,iq)) CALL WriteField_u('zm',zm(:,:,iq)) #endif SELECT CASE(tracers(iq)%iadv) CASE(0); CYCLE CASE(10); call vly_loc(zq,pente_max,zm,mv, iq) CASE(14); call vlyqs_loc(zq,pente_max,zm,mv,qsat,iq) CASE DEFAULT CALL abort_gcm("vlspltgen_p","schema non parallelise",1) END SELECT enddo !do iq=1,nqtot call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 493') do iq=1,nqtot IF(tracers(iq)%parent /= 'air') CYCLE ! !write(*,*) 'vlspltgen 477: iq=',iq #ifdef DEBUG_IO CALL WriteField_u('zq',zq(:,:,iq)) CALL WriteField_u('zm',zm(:,:,iq)) #endif SELECT CASE(tracers(iq)%iadv) CASE(0); CYCLE CASE(10); call vlx_loc(zq,pente_max,zm,mu, & ij_begin,ij_end,iq) CASE(14); call vlxqs_loc(zq,pente_max,zm,mu, & qsat, ij_begin,ij_end,iq) CASE DEFAULT CALL abort_gcm("vlspltgen_p","schema non parallelise",1) END SELECT enddo !do iq=1,nqtot ! !write(*,*) 'vlspltgen 550: apres derniere serie de call vlx' call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 521') ijb=ij_begin ije=ij_end ! !write(*,*) 'vlspltgen_loc 557' !$OMP BARRIER ! !write(*,*) 'vlspltgen_loc 559' DO iq=1,nqtot ! !write(*,*) 'vlspltgen_loc 561, iq=',iq #ifdef DEBUG_IO CALL WriteField_u('zq',zq(:,:,iq)) CALL WriteField_u('zm',zm(:,:,iq)) #endif !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm DO ij=ijb,ije ! print *,'zq-->',ij,l,iq,zq(ij,l,iq) ! print *,'q-->',ij,l,iq,q(ij,l,iq) q(ij,l,iq)=zq(ij,l,iq) ENDDO ENDDO !$OMP END DO NOWAIT ! !write(*,*) 'vlspltgen_loc 575' !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,llm DO ij=ijb,ije-iip1+1,iip1 q(ij+iim,l,iq)=q(ij,l,iq) ENDDO ENDDO !$OMP END DO NOWAIT ! !write(*,*) 'vlspltgen_loc 583' ENDDO !DO iq=1,nqtot call check_isotopes(q,ij_begin,ij_end,'vlspltgen_loc 557') !$OMP BARRIER !c$OMP MASTER ! call WaitSendRequest(MyRequest1) ! call WaitSendRequest(MyRequest2) !c$OMP END MASTER !c$OMP BARRIER ! !write(*,*) 'vlspltgen 597: sortie' RETURN END SUBROUTINE vlspltgen_loc