MODULE lmdz_filtreg_p USE lmdz_filtreg, ONLY: matrinvn, matrinvs, matriceun, matriceus, matricevn, matricevs IMPLICIT NONE; PRIVATE PUBLIC filtreg_p CONTAINS SUBROUTINE filtreg_p(champ, jjb, jje, ibeg, iend, nlat, nbniv, & ifiltre, iaire, griscal, iter) USE parallel_lmdz, ONLY: OMP_CHUNK USE lmdz_filtre_fft_loc, ONLY: use_filtre_fft, filtre_u_fft, & filtre_v_fft, filtre_inv_fft USE lmdz_timer_filtre, ONLY: init_timer, start_timer, stop_timer USE lmdz_coefils, ONLY: jfiltnu, jfiltnv, jfiltsu, jfiltsv, sddu, sddv, unsddu, unsddv, modfrstv, modfrstu !======================================================================= ! ! Auteur: P. Le Van 07/10/97 ! ------ ! ! Objet: filtre matriciel longitudinal ,avec les matrices precalculees ! pour l'operateur Filtre . ! ------ ! ! Arguments: ! ---------- ! ! ! ibeg..iend lattitude a filtrer ! nlat nombre de latitudes du champ ! nbniv nombre de niveaux verticaux a filtrer ! champ(iip1,nblat,nbniv) en entree : champ a filtrer ! en sortie : champ filtre ! ifiltre +1 Transformee directe ! -1 Transformee inverse ! +2 Filtre directe ! -2 Filtre inverse ! ! iaire 1 si champ intensif ! 2 si champ extensif (pondere par les aires) ! ! iter 1 filtre simple ! !======================================================================= ! ! ! Variable Intensive ! ifiltre = 1 filtre directe ! ifiltre =-1 filtre inverse ! ! Variable Extensive ! ifiltre = 2 filtre directe ! ifiltre =-2 filtre inverse ! ! INCLUDE "dimensions.h" INCLUDE "paramet.h" ! INTEGER, INTENT(IN) :: jjb, jje, ibeg, iend, nlat, nbniv, ifiltre, iter INTEGER, INTENT(IN) :: iaire LOGICAL, INTENT(IN) :: griscal REAL, INTENT(INOUT) :: champ(iip1, jjb:jje, nbniv) INTEGER :: i, j, l, k INTEGER :: iim2, immjm INTEGER :: jdfil1, jdfil2, jffil1, jffil2, jdfil, jffil INTEGER :: hemisph REAL :: champ_fft(iip1, jjb:jje, nbniv) ! REAL :: champ_in(iip1,jjb:jje,nbniv) LOGICAL, SAVE :: first = .TRUE. !$OMP THREADPRIVATE(first) REAL, DIMENSION(iip1, jjb:jje, nbniv) :: champ_loc INTEGER :: ll_nb, nbniv_loc REAL, SAVE :: sdd12(iim, 4) !$OMP THREADPRIVATE(sdd12) INTEGER, PARAMETER :: type_sddu = 1 INTEGER, PARAMETER :: type_sddv = 2 INTEGER, PARAMETER :: type_unsddu = 3 INTEGER, PARAMETER :: type_unsddv = 4 INTEGER :: sdd1_type, sdd2_type CHARACTER (LEN = 132) :: abort_message IF (first) THEN sdd12(1:iim, type_sddu) = sddu(1:iim) sdd12(1:iim, type_sddv) = sddv(1:iim) sdd12(1:iim, type_unsddu) = unsddu(1:iim) sdd12(1:iim, type_unsddv) = unsddv(1:iim) CALL Init_timer first = .FALSE. ENDIF !$OMP MASTER CALL start_timer !$OMP END MASTER !-------------------------------------------------------c IF(ifiltre==1.OR.ifiltre==-1) & CALL abort_gcm("lmdz_filtreg_p", 'Pas de transformee& &simple dans cette version', 1) IF(iter== 2) THEN PRINT *, ' Pas d iteration du filtre dans cette version !'& &, ' Utiliser old_filtreg et repasser !' CALL abort_gcm("lmdz_filtreg_p", "stopped", 1) ENDIF IF(ifiltre== -2 .AND..NOT.griscal) THEN PRINT *, ' Cette routine ne calcule le filtre inverse que ' & , ' sur la grille des scalaires !' CALL abort_gcm("lmdz_filtreg_p", "stopped", 1) ENDIF IF(ifiltre/=2 .AND.ifiltre/= - 2) THEN PRINT *, ' Probleme dans filtreg car ifiltre NE 2 et NE -2' & , ' corriger et repasser !' CALL abort_gcm("lmdz_filtreg_p", "stopped", 1) ENDIF ! iim2 = iim * iim immjm = iim * jjm ! ! IF(griscal) THEN IF(nlat /= jjp1) THEN CALL abort_gcm("lmdz_filtreg_p", " nlat. NE. jjp1", 1) ELSE ! IF(iaire==1) THEN sdd1_type = type_sddv sdd2_type = type_unsddv ELSE sdd1_type = type_unsddv sdd2_type = type_sddv ENDIF ! jdfil1 = 2 jffil1 = jfiltnu jdfil2 = jfiltsu jffil2 = jjm ENDIF ELSE IF(nlat/=jjm) THEN CALL abort_gcm("lmdz_filtreg_p", " nlat. NE. jjm", 1) ELSE ! IF(iaire==1) THEN sdd1_type = type_sddu sdd2_type = type_unsddu ELSE sdd1_type = type_unsddu sdd2_type = type_sddu ENDIF ! jdfil1 = 1 jffil1 = jfiltnv jdfil2 = jfiltsv jffil2 = jjm ENDIF ENDIF ! DO hemisph = 1, 2 ! IF (hemisph==1) THEN !ym jdfil = max(jdfil1, ibeg) jffil = min(jffil1, iend) ELSE !ym jdfil = max(jdfil2, ibeg) jffil = min(jffil2, iend) ENDIF !ccccccccccccccccccccccccccccccccccccccccccc ! Utilisation du filtre classique !ccccccccccccccccccccccccccccccccccccccccccc IF (.NOT. use_filtre_fft) THEN !---------------------------------! ! Agregation des niveau verticaux ! ! uniquement necessaire pour une ! ! execution OpenMP ! !---------------------------------! ll_nb = 0 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, nbniv ll_nb = ll_nb + 1 DO j = jdfil, jffil DO i = 1, iim champ_loc(i, j, ll_nb) = & champ(i, j, l) * sdd12(i, sdd1_type) ENDDO ENDDO ENDDO !$OMP END DO NOWAIT nbniv_loc = ll_nb IF(hemisph==1) THEN IF(ifiltre==-2) THEN DO j = jdfil, jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, & matrinvn(1,1,j), iim, & champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, & champ_fft(1,j,1), iip1*(jje-jjb+1)) #else champ_fft(1:iim, j, 1:nbniv_loc) = & matmul(matrinvn(1:iim, 1:iim, j), & champ_loc(1:iim, j, 1:nbniv_loc)) #endif ENDDO ELSE IF (griscal) THEN DO j = jdfil, jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, & matriceun(1,1,j), iim, & champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, & champ_fft(1,j,1), iip1*(jje-jjb+1)) #else champ_fft(1:iim, j, 1:nbniv_loc) = & matmul(matriceun(1:iim, 1:iim, j), & champ_loc(1:iim, j, 1:nbniv_loc)) #endif ENDDO ELSE DO j = jdfil, jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, & matricevn(1,1,j), iim, & champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, & champ_fft(1,j,1), iip1*(jje-jjb+1)) #else champ_fft(1:iim, j, 1:nbniv_loc) = & matmul(matricevn(1:iim, 1:iim, j), & champ_loc(1:iim, j, 1:nbniv_loc)) #endif ENDDO ENDIF ELSE IF(ifiltre==-2) THEN DO j = jdfil, jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, & matrinvs(1,1,j-jfiltsu+1), iim, & champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, & champ_fft(1,j,1), iip1*(jje-jjb+1)) #else champ_fft(1:iim, j, 1:nbniv_loc) = & matmul(matrinvs(1:iim, 1:iim, j - jfiltsu + 1), & champ_loc(1:iim, j, 1:nbniv_loc)) #endif ENDDO ELSE IF (griscal) THEN DO j = jdfil, jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, & matriceus(1,1,j-jfiltsu+1), iim, & champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, & champ_fft(1,j,1), iip1*(jje-jjb+1)) #else champ_fft(1:iim, j, 1:nbniv_loc) = & matmul(matriceus(1:iim, 1:iim, j - jfiltsu + 1), & champ_loc(1:iim, j, 1:nbniv_loc)) #endif ENDDO ELSE DO j = jdfil, jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, & matricevs(1,1,j-jfiltsv+1), iim, & champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, & champ_fft(1,j,1), iip1*(jje-jjb+1)) #else champ_fft(1:iim, j, 1:nbniv_loc) = & matmul(matricevs(1:iim, 1:iim, j - jfiltsv + 1), & champ_loc(1:iim, j, 1:nbniv_loc)) #endif ENDDO ENDIF ENDIF ! c IF(ifiltre==2) THEN !-------------------------------------! ! Dés-agregation des niveau verticaux ! ! uniquement necessaire pour une ! ! execution OpenMP ! !-------------------------------------! ll_nb = 0 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, nbniv ll_nb = ll_nb + 1 DO j = jdfil, jffil DO i = 1, iim champ(i, j, l) = (champ_loc(i, j, ll_nb) & + champ_fft(i, j, ll_nb)) & * sdd12(i, sdd2_type) ENDDO ENDDO ENDDO !$OMP END DO NOWAIT ELSE ll_nb = 0 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, nbniv ll_nb = ll_nb + 1 DO j = jdfil, jffil DO i = 1, iim champ(i, j, l) = (champ_loc(i, j, ll_nb) & - champ_fft(i, j, ll_nb)) & * sdd12(i, sdd2_type) ENDDO ENDDO ENDDO !$OMP END DO NOWAIT ENDIF !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, nbniv DO j = jdfil, jffil ! add redundant longitude champ(iip1, j, l) = champ(1, j, l) ENDDO ENDDO !$OMP END DO NOWAIT !cccccccccccccccccccccccccccccccccccccccccccc ! Utilisation du filtre FFT !cccccccccccccccccccccccccccccccccccccccccccc ELSE !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, nbniv DO j = jdfil, jffil DO i = 1, iim champ(i, j, l) = champ(i, j, l) * sdd12(i, sdd1_type) champ_fft(i, j, l) = champ(i, j, l) ENDDO ENDDO ENDDO !$OMP END DO NOWAIT IF (jdfil<=jffil) THEN IF(ifiltre == -2) THEN CALL Filtre_inv_fft(champ_fft, jjb, jje, jdfil, jffil, nbniv) ELSE IF (griscal) THEN CALL Filtre_u_fft(champ_fft, jjb, jje, jdfil, jffil, nbniv) ELSE CALL Filtre_v_fft(champ_fft, jjb, jje, jdfil, jffil, nbniv) ENDIF ENDIF IF(ifiltre== 2) THEN !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, nbniv DO j = jdfil, jffil DO i = 1, iim champ(i, j, l) = (champ(i, j, l) + champ_fft(i, j, l)) & * sdd12(i, sdd2_type) ENDDO ENDDO ENDDO !$OMP END DO NOWAIT ELSE !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, nbniv DO j = jdfil, jffil DO i = 1, iim champ(i, j, l) = (champ(i, j, l) - champ_fft(i, j, l)) & * sdd12(i, sdd2_type) ENDDO ENDDO ENDDO !$OMP END DO NOWAIT ENDIF ! !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, nbniv DO j = jdfil, jffil ! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l ) ! ! add redundant longitude champ(iip1, j, l) = champ(1, j, l) ENDDO ENDDO !$OMP END DO NOWAIT ENDIF ! Fin de la zone de filtrage ENDDO ! DO j=1,nlat ! PRINT *,"check FFT ----> Delta(",j,")=", ! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)), ! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:)) ! ENDDO ! PRINT *,"check FFT ----> Delta(",j,")=", ! & sum(champ-champ_fft)/sum(champ) ! !$OMP MASTER CALL stop_timer !$OMP END MASTER END SUBROUTINE filtreg_p END MODULE lmdz_filtreg_p