From 5ea698e1cbbbeb6cd3833bd884f8d65e899e5d83 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Wed, 17 Dec 2025 14:43:36 +0000 Subject: [PATCH 1/7] Migrate um13.9_chemistry_optim over from FCM --- .../core/chemistry/asad/asad_cdrive.F90 | 39 +- src/science/core/chemistry/asad/asad_ftoy.F90 | 50 +-- src/science/core/chemistry/asad/asad_prls.F90 | 6 +- .../core/chemistry/asad/asad_sparse_vars.F90 | 345 ++++++++---------- .../core/chemistry/asad/asad_spimpmjp.F90 | 5 +- .../core/chemistry/asad/asad_steady.F90 | 246 ++++++------- .../core/chemistry/ukca_chemistry_ctl.F90 | 8 +- .../chemistry/ukca_chemistry_ctl_col_mod.F90 | 7 +- .../chemistry/ukca_chemistry_ctl_full_mod.F90 | 5 +- 9 files changed, 283 insertions(+), 428 deletions(-) diff --git a/src/science/core/chemistry/asad/asad_cdrive.F90 b/src/science/core/chemistry/asad/asad_cdrive.F90 index 17b6750..b3f0629 100644 --- a/src/science/core/chemistry/asad/asad_cdrive.F90 +++ b/src/science/core/chemistry/asad/asad_cdrive.F90 @@ -33,7 +33,6 @@ ! This file belongs in section: UKCA ! ! Arguments: -! cdot - Tracer tendencies due to chemistry. ! ftr - Tracer concentrations. ! pp - Pressure (Nm-2). ! pt - Temperature (K). @@ -90,11 +89,11 @@ MODULE asad_cdrive_mod CONTAINS -SUBROUTINE asad_cdrive(cdot, ftr, pp, pt, pq, co2_1d, cld_f, cld_l, & +SUBROUTINE asad_cdrive(ftr, pp, pt, pq, co2_1d, cld_f, cld_l, & ix, jy, nlev, dryrt, wetrt, rc_het, prt, & n_points, have_nat, stratflag, H_plus_1d_arr) -USE asad_mod, ONLY: cdt, ctype, fdot, f, & +USE asad_mod, ONLY: ctype, fdot, f, & jpspec, jpcspf, jppj, & jpdd, jpdw, jpif, jsubs, & linfam, lvmr, & @@ -155,8 +154,6 @@ SUBROUTINE asad_cdrive(cdot, ftr, pp, pt, pq, co2_1d, cld_f, cld_l, & ! rather than number of tracers REAL, INTENT(IN OUT) :: ftr(n_points,jpcspf) ! Tracer concs -REAL, INTENT(OUT) :: cdot(n_points,jpcspf) ! Tracer tendencies - ! Local variables INTEGER :: errcode ! Variable passed to ereport @@ -187,24 +184,16 @@ SUBROUTINE asad_cdrive(cdot, ftr, pp, pt, pq, co2_1d, cld_f, cld_l, & ! 1. Initialise variables and arrays -! 1.1 Clear tendencies to avoid contributions from levels -! on which no chemistry is performed - IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) -DO jtr = 1, jpcspf - DO jl = 1, n_points - cdot(jl,jtr) = 0.0 - END DO -END DO -! 1.2 Copy pressure and temperature to asad_mod +! 1.1 Copy pressure and temperature to asad_mod DO jl = 1, n_points p(jl) = pp(jl) t(jl) = pt(jl) END DO -! 1.2.1 Copy water vapor and co2 to asad_mod +! 1.1.1 Copy water vapor and co2 to asad_mod DO jl = 1, n_points wp(jl) = pq(jl) @@ -389,8 +378,7 @@ SUBROUTINE asad_cdrive(cdot, ftr, pp, pt, pq, co2_1d, cld_f, cld_l, & itr = madvtr(js) iodd = nodd(js) DO jl = 1, n_points - IF ( linfam(jl,itr) ) f(jl,ifam) = & - f(jl,ifam) - iodd*f(jl,itr) + IF ( linfam(jl,itr) ) f(jl,ifam) = f(jl,ifam) - iodd*f(jl,itr) END DO END IF END DO @@ -398,17 +386,9 @@ SUBROUTINE asad_cdrive(cdot, ftr, pp, pt, pq, co2_1d, cld_f, cld_l, & ! 7.2 Returned values of concentration and chemical tendency DO jtr = 1, jpcspf - IF ( method /= 0 ) THEN - DO jl = 1, n_points - cdot(jl,jtr) = ( f(jl,jtr)-ftr(jl,jtr)) / (cdt*ncsteps) - ftr(jl,jtr) = f(jl,jtr) - END DO - ELSE - DO jl = 1, n_points - cdot(jl,jtr) = fdot(jl,jtr) - ftr(jl,jtr) = f(jl,jtr) - END DO - END IF + DO jl = 1, n_points + ftr(jl,jtr) = f(jl,jtr) + END DO END DO @@ -418,8 +398,7 @@ SUBROUTINE asad_cdrive(cdot, ftr, pp, pt, pq, co2_1d, cld_f, cld_l, & IF ( lvmr ) THEN DO jtr = 1, jpcspf DO jl = 1, n_points - ftr(jl,jtr) = ftr(jl,jtr) / tnd(jl) - cdot(jl,jtr) = cdot(jl,jtr) / tnd(jl) + ftr(jl,jtr) = ftr(jl,jtr) / tnd(jl) END DO END DO END IF diff --git a/src/science/core/chemistry/asad/asad_ftoy.F90 b/src/science/core/chemistry/asad/asad_ftoy.F90 index b451263..43ebaca 100644 --- a/src/science/core/chemistry/asad/asad_ftoy.F90 +++ b/src/science/core/chemistry/asad/asad_ftoy.F90 @@ -182,7 +182,6 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) INTEGER :: js ! Loop variable INTEGER :: ifam INTEGER :: imaj -INTEGER :: nl INTEGER :: istart INTEGER :: iend INTEGER :: iodd @@ -190,11 +189,8 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) INTEGER :: icode ! Error code INTEGER :: iro2 ! Counter for RO2 species - INTEGER, SAVE :: istmin INTEGER, SAVE :: ift -! INTEGER :: ilstmin(jpspec) ! Now in ASAD_MOD -! INTEGER :: ilft(jpspec) ! Now in ASAD_MOD CHARACTER (LEN=errormessagelength) :: cmessage ! Error message INTEGER :: errcode ! Variable passed to ereport @@ -222,8 +218,8 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) ! 1. Initialise. ! ----------- -IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName, & - zhook_in,zhook_handle) +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) + ! OMP CRITICAL will only allow one thread through this code at a time, ! while the other threads are held until completion. !$OMP CRITICAL (asad_ftoy_init) @@ -291,7 +287,7 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) DO j = 1, jpro2 ! Get index location of each RO2 species and sum - iro2 = nlfro2(j) + iro2 = nlfro2(j) fro2(1:n_points) = fro2(1:n_points) + f(1:n_points, iro2) END DO ! End iteration over RO2 species @@ -311,23 +307,19 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) END IF ! Calculate the new concentration of RO2 - y(1:n_points, jsro2) = fro2_old(1:n_points) + & - damp_ro2*(fro2(1:n_points) - fro2_old(1:n_points)) + y(1:n_points, jsro2) = fro2_old(1:n_points) + damp_ro2*(fro2(1:n_points) - fro2_old(1:n_points)) END IF ! 1.2 If there are no family or steady species then exit IF ( nstst == 0) THEN - IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName, & - zhook_out,zhook_handle) + IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN END IF ! 1.3 Initialise local variables and do sanity checks. -nl = n_points - DO j = 1, nstst js = nlstst(j) DO jl = 1, n_points @@ -354,8 +346,7 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) IF (iter == 0) THEN CALL asad_fyfixr(n_points) - IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName, & - zhook_out,zhook_handle) + IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN END IF @@ -373,10 +364,7 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) DO j = 1, nstst js = nlstst(j) DO jl = 1, n_points - IF (ABS(y(jl,js)) < peps .AND. & - prod(jl,js) > peps ) THEN - y(jl,js) = peps - END IF + IF (ABS(y(jl,js)) < peps .AND. prod(jl,js) > peps ) y(jl,js) = peps END DO END DO END IF @@ -387,7 +375,7 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) ! Regardless of whether the user has them turned on or not. ! See method above. - CALL asad_prls( nl, istmin, ilstmin, gdepem ) + CALL asad_prls( n_points, istmin, ilstmin, gdepem ) ! 4.2 Initialise ratios @@ -446,8 +434,7 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) DO jl = 1, n_points IF ( y(jl,imaj) > peps ) THEN ratio(jl,js) = y(jl,js) / y(jl,imaj) - ratio(jl,imaj) = ratio(jl,imaj) + & - iodd * ratio(jl,js) + ratio(jl,imaj) = ratio(jl,imaj) + iodd * ratio(jl,js) ELSE ratio(jl,js) = 0.0 END IF @@ -471,9 +458,8 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) ELSE sl = 0.0 END IF - linfam(jl,itr) = sl > zthresh - IF ( linfam(jl,itr) ) & - f(jl,ifam) = f(jl,ifam) + iodd*f(jl,itr) + linfam(jl,itr) = sl > zthresh + IF ( linfam(jl,itr) ) f(jl,ifam) = f(jl,ifam) + iodd*f(jl,itr) END DO END IF @@ -481,8 +467,7 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) IF ( linfam(jl,itr) ) THEN IF ( y(jl,imaj) > peps ) THEN ratio(jl,js) = y(jl,js) / y(jl,imaj) - ratio(jl,imaj) = ratio(jl,imaj) + & - iodd * ratio(jl,js) + ratio(jl,imaj) = ratio(jl,imaj) + iodd * ratio(jl,js) ELSE ratio(jl,js) = 0.0 END IF @@ -520,8 +505,7 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) ELSE itr = madvtr(js) DO jl = 1, n_points - IF ( linfam(jl,itr) ) y(jl,js) = & - y(jl,imaj) * ratio(jl,js) + IF ( linfam(jl,itr) ) y(jl,js) = y(jl,imaj) * ratio(jl,js) END DO END IF END DO @@ -538,8 +522,7 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) IF ( .NOT. gconv ) EXIT END DO IF (gconv) THEN - IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName, & - zhook_out,zhook_handle) + IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN END IF @@ -552,11 +535,10 @@ SUBROUTINE asad_ftoy(ofirst,iter, num_iter, n_points, ix, jy, nlev) END DO ! End of iterations - ! 9. Convergence achieved or max. iterations reached. -IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName, & - zhook_out,zhook_handle) +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN END SUBROUTINE asad_ftoy + END MODULE asad_ftoy_mod diff --git a/src/science/core/chemistry/asad/asad_prls.F90 b/src/science/core/chemistry/asad/asad_prls.F90 index 6d7613a..315de77 100644 --- a/src/science/core/chemistry/asad/asad_prls.F90 +++ b/src/science/core/chemistry/asad/asad_prls.F90 @@ -219,8 +219,7 @@ SUBROUTINE asad_prls( kl, knspec, kspec, ldepem ) i3 = nprdx3(3,j3,js) DO jl = 1, kl - pd(jl,js) = pd(jl,js) & - + prk(jl,i1) + prk(jl,i2) + prk(jl,i3) + pd(jl,js) = pd(jl,js) + prk(jl,i1) + prk(jl,i2) + prk(jl,i3) END DO END DO @@ -271,8 +270,7 @@ SUBROUTINE asad_prls( kl, knspec, kspec, ldepem ) DO j = istart, iend js = nldepx(j) DO jl = 1, kl - slos(jl,js) = slos(jl,js) + ( dpd(jl,js) + dpw(jl,js) ) & - * y(jl,js) + slos(jl,js) = slos(jl,js) + ( dpd(jl,js) + dpw(jl,js) ) * y(jl,js) END DO END DO diff --git a/src/science/core/chemistry/asad/asad_sparse_vars.F90 b/src/science/core/chemistry/asad/asad_sparse_vars.F90 index 01aacc9..8bc69df 100644 --- a/src/science/core/chemistry/asad/asad_sparse_vars.F90 +++ b/src/science/core/chemistry/asad/asad_sparse_vars.F90 @@ -142,16 +142,14 @@ SUBROUTINE setup_spfuljac() INTEGER :: kr INTEGER :: ikr INTEGER :: krj -INTEGER :: ij1 +INTEGER :: ij INTEGER :: itemp1 -INTEGER :: ij(jpmsp) INTEGER :: activity(jpcspf) INTEGER, ALLOCATABLE :: permute(:,:) ! permutation matrix INTEGER, ALLOCATABLE :: map(:,:) ! a matrix of 1's indicating where the ! dense Jacobian A(i,j) has nonzero ! entries -LOGICAL :: lj(jpmsp) CHARACTER(LEN=errormessagelength) :: cmessage @@ -169,22 +167,20 @@ SUBROUTINE setup_spfuljac() ! ------------------------------------------------------------------------------ IF (.NOT. ALLOCATED(map)) ALLOCATE(map(jpcspf, jpcspf)) + map(:,:) = 0 -DO i=1,jpcspf +DO i = 1, jpcspf map(i,i) = 1 END DO DO jc = 1, ntrf itrcr = nltrf(jc) - DO j3 = 1,nmzjac(itrcr) + DO j3 = 1, nmzjac(itrcr) irj = nzjac1(j3,itrcr) ! - DO jn=1,jpmsp - ij(jn)=njcoth(irj,jn) - lj(jn)= (ij(jn) /= 0) - END DO - DO jn=1,jpmsp - IF (lj(jn)) map(ij(jn),itrcr) = 1 + DO jn = 1, jpmsp + ij = njcoth(irj,jn) + IF (ij /= 0) map(ij,itrcr) = 1 END DO ! IF (npdfr(irj,1) /= 0) THEN @@ -204,19 +200,11 @@ SUBROUTINE setup_spfuljac() ! required for the O3 loss rate. ! DO jc = 1, nstst - DO j3 = 1,nmsjac(jc) + DO j3 = 1, nmsjac(jc) irj = nsjac1(j3,jc) - DO jn=1,jpmsp - ij(jn)=njcoth(irj,jn) - lj(jn)= (ij(jn) /= 0) - END DO - DO jn=1,jpmsp - IF (lj(jn)) THEN - map(ij(jn),ntro3 ) = 1 - ! map(ij(jn),ntroh ) = 1 - ! map(ij(jn),ntrho2) = 1 - ! map(ij(jn),ntrno ) = 1 - END IF + DO jn = 1, jpmsp + ij = njcoth(irj,jn) + IF (ij /= 0) map(ij,ntro3) = 1 END DO END DO END DO @@ -225,10 +213,10 @@ SUBROUTINE setup_spfuljac() ! Add deposition terms to Jacobian diagonal. ! ----------------------------------------------------------------- ! -IF ( (ndepw /= 0) .OR. (ndepd /= 0) ) THEN +IF ((ndepw /= 0) .OR. (ndepd /= 0) ) THEN DO js = 1, jpspec itrd = madvtr(js) - IF ( itrd /= 0 ) map(itrd,itrd) = 1 + IF (itrd /= 0) map(itrd,itrd) = 1 END DO END IF @@ -239,10 +227,10 @@ SUBROUTINE setup_spfuljac() CALL umPrint(umMessage,src='asad_sparse_vars') END IF -isx=0 -nonzero_map(:,:)=0 -DO i=1,jpcspf - DO j=1,jpcspf +isx = 0 +nonzero_map(:,:) = 0 +DO i = 1, jpcspf + DO j = 1, jpcspf ! calculate forward and backward pointers for sparse representation IF (map(i,j) == 1) THEN isx = isx + 1 @@ -253,10 +241,10 @@ SUBROUTINE setup_spfuljac() END DO base_tracer(:) = 0 -DO i=1,jpcspf - DO j=1,jpcspf +DO i = 1, jpcspf + DO j = 1, jpcspf i1 = nonzero_map(j,i) - IF (i1 > 0) base_tracer(i1)=i + IF (i1 > 0) base_tracer(i1) = i END DO END DO @@ -264,14 +252,14 @@ SUBROUTINE setup_spfuljac() ! Section 2: Create nonzero_map_unordered (index array for PAP') ! ------------------------------------------------------------------------------ -! Reorder species to minimize fill-in -DO i=1,jpcspf +! Reorder species by their reactivity to minimize fill-in +DO i = 1, jpcspf activity(i) = SUM(map(:,i)) + SUM(map(i,:)) reorder(i) = i END DO -DO i=1,jpcspf-1 - DO j=i+1,jpcspf +DO i = 1, jpcspf-1 + DO j = i+1, jpcspf IF (activity(i) > activity(j)) THEN ! exchange i and j tracers if i is more active than j. itemp1 = reorder(i) @@ -284,7 +272,7 @@ SUBROUTINE setup_spfuljac() END DO END DO -DO i=1,jpcspf +DO i = 1, jpcspf IF (mype == 0 .AND. printstatus >= prstatus_oper) THEN WRITE(umMessage,'(A,I3,1X,I3,1X,A10,1X,I3)') 'IN ASAD SETUP_SPFULJAC: ', & i, reorder(i), specf(reorder(i)), activity(i) @@ -296,7 +284,7 @@ SUBROUTINE setup_spfuljac() IF (.NOT. ALLOCATED(permute)) ALLOCATE(permute(jpcspf,jpcspf)) permute(:,:) = 0 -DO i=1,jpcspf +DO i = 1, jpcspf permute(i,reorder(i)) = 1 END DO @@ -304,10 +292,10 @@ SUBROUTINE setup_spfuljac() ! spfj(nonzero_map_unordered(i,j)) = P*A*P'(i,j) ! where P is the permutation matrix and A is the dense Jacobian matrix ! with spfj the compressed storage (sparse) array for matrix A. -nonzero_map_unordered = MATMUL(MATMUL(permute,nonzero_map),TRANSPOSE(permute)) +nonzero_map_unordered = MATMUL(MATMUL(permute, nonzero_map), TRANSPOSE(permute)) IF (ALLOCATED(permute)) DEALLOCATE(permute) -IF (ALLOCATED(map)) DEALLOCATE(map) +IF (ALLOCATED(map)) DEALLOCATE(map) ! ------------------------------------------------------------------------------ ! Section 3: Check the number of nonzero elements in LU factorization @@ -317,17 +305,17 @@ SUBROUTINE setup_spfuljac() ! array PAP' and check that it is less than spfjsize_max. total1 = total modified_map(:,:) = nonzero_map_unordered(:,:) -DO kr=1,jpcspf - DO i=kr+1,jpcspf +DO kr = 1, jpcspf + DO i = kr+1, jpcspf ikr = modified_map(i,kr) IF (ikr > 0) THEN - DO j=kr+1,jpcspf + DO j = kr+1, jpcspf krj = modified_map(kr,j) IF (krj > 0) THEN - ij1 = modified_map(i,j) + ij = modified_map(i,j) ! Distinguish whether matrix element is zero or not. If not, proceed ! as in dense case. If it is, create new matrix element. - IF (ij1 <= 0) THEN + IF (ij <= 0) THEN total1 = total1 + 1 modified_map(i,j) = total1 END IF @@ -375,54 +363,44 @@ SUBROUTINE setup_spfuljac() DO jc = 1, ntrf itrcr = nltrf(jc) ! - DO j3 = 1,nmzjac(itrcr) + DO j3 = 1, nmzjac(itrcr) irj = nzjac1(j3,itrcr) ! - DO jn=1,jpmsp - ij(jn)=njcoth(irj,jn) - lj(jn)= (ij(jn) /= 0) - END DO - DO jn=1,2 - IF (lj(jn)) THEN - i = nonzero_map(ij(jn),itrcr) - nnegterms(i) = nnegterms(i) + 1 - IF (nnegterms(i) > maxterms) THEN - errcodes(i,1) = 1 - ELSE IF (.NOT. ANY(errcodes(i,:) /= 0)) THEN - negterms(i,nnegterms(i)) = irj - END IF - END IF - END DO - IF (nfrpx(irj) == 0) THEN - DO jn=3,jpmsp - IF (lj(jn)) THEN - i = nonzero_map(ij(jn),itrcr) - nposterms(i) = nposterms(i) + 1 - IF (nposterms(i) > maxterms) THEN - errcodes(i,1) = 2 - ELSE IF (.NOT. ANY(errcodes(i,:) /= 0)) THEN - posterms(i,nposterms(i)) = irj - END IF - END IF - END DO - ELSE - DO jn=3,jpmsp - IF (lj(jn)) THEN - i = nonzero_map(ij(jn),itrcr) - nfracterms(i) = nfracterms(i) + 1 - IF (nfracterms(i) > maxfterms) THEN - errcodes(i,1) = 3 + DO jn = 1, jpmsp + ij = njcoth(irj,jn) + IF (ij /= 0) THEN + i = nonzero_map(ij,itrcr) + IF (jn < 3) THEN + nnegterms(i) = nnegterms(i) + 1 + IF (nnegterms(i) > maxterms) THEN + errcodes(i,1) = 1 + ELSE IF (ALL(errcodes(i,:) == 0)) THEN + negterms(i,nnegterms(i)) = irj END IF - IF (nfrpx(irj)+jn-3 > jpfrpx) THEN - errcodes(i,2) = i - errcodes(i,3) = jn - ELSE IF (.NOT. ANY(errcodes(i,:) /= 0)) THEN - fracterms(i,nfracterms(i)) = irj - ffrac(i,nfracterms(i)) = frpx(nfrpx(irj)+jn-3) + ELSE + IF (nfrpx(irj) == 0) THEN + nposterms(i) = nposterms(i) + 1 + IF (nposterms(i) > maxterms) THEN + errcodes(i,1) = 2 + ELSE IF (ALL(errcodes(i,:) == 0)) THEN + posterms(i,nposterms(i)) = irj + END IF + ELSE + nfracterms(i) = nfracterms(i) + 1 + IF (nfracterms(i) > maxfterms) THEN + errcodes(i,1) = 3 + END IF + IF (nfrpx(irj)+jn-3 > jpfrpx) THEN + errcodes(i,2) = i + errcodes(i,3) = jn + ELSE IF (ALL(errcodes(i,:) == 0)) THEN + fracterms(i,nfracterms(i)) = irj + ffrac(i,nfracterms(i)) = frpx(nfrpx(irj)+jn-3) + END IF END IF END IF - END DO - END IF + END IF + END DO ! IF (npdfr(irj,1) /= 0) THEN i1 = npdfr(irj,1) @@ -471,7 +449,7 @@ END SUBROUTINE setup_spfuljac !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -SUBROUTINE spfuljac(n_points,cdt,min_pivot,nonzero_map,spfj) +SUBROUTINE spfuljac(n_points, cdt, min_pivot, nonzero_map, spfj) ! ! Routine to calculate the Jacobian in sparse format ! @@ -507,9 +485,7 @@ SUBROUTINE spfuljac(n_points,cdt,min_pivot,nonzero_map,spfj) INTEGER :: jn INTEGER :: js INTEGER :: jl -INTEGER :: ij(jpmsp) - -LOGICAL :: lj(jpmsp) +INTEGER :: ij INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 @@ -521,21 +497,20 @@ SUBROUTINE spfuljac(n_points,cdt,min_pivot,nonzero_map,spfj) ! At the bottom of this routine we divide by f, so ensure that ! f is not too small. -DO jl = 1, jpcspf - DO i = 1, n_points - IF (f(i, jl) < min_pivot) THEN - f(i, jl) = min_pivot +DO i = 1, jpcspf + DO jl = 1, n_points + IF (f(jl, i) < min_pivot) THEN + f(jl, i) = min_pivot END IF END DO END DO -deltt=1.0/cdt - -spfj(:,:) = 0.0 +deltt = 1.0 / cdt ! Calculate diagonal element of Jacobian -DO i=1,jpcspf - spfj(:,nonzero_map(i,i))=-deltt*f(:,i) +spfj(:,:) = 0.0 +DO i = 1, jpcspf + spfj(:,nonzero_map(i,i)) = -deltt*f(:,i) END DO ! ! @@ -545,14 +520,14 @@ SUBROUTINE spfuljac(n_points,cdt,min_pivot,nonzero_map,spfj) ! ! Sum up positive non-fractional, negative, and fractional terms -DO p=1,total - DO i=nposterms(p),1,-1 +DO p = 1, total + DO i = nposterms(p), 1, -1 spfj(:,p) = spfj(:,p) + prk(:,posterms(p,i)) END DO - DO i=nnegterms(p),1,-1 + DO i = nnegterms(p), 1, -1 spfj(:,p) = spfj(:,p) - prk(:,negterms(p,i)) END DO - DO i=nfracterms(p),1,-1 + DO i = nfracterms(p), 1, -1 spfj(:,p) = spfj(:,p) + ffrac(p,i)*prk(:,fracterms(p,i)) END DO END DO @@ -562,24 +537,17 @@ SUBROUTINE spfuljac(n_points,cdt,min_pivot,nonzero_map,spfj) ! required for the O3 loss rate. ! DO jc = 1, nstst - DO j3 = 1,nmsjac(jc) + DO j3 = 1, nmsjac(jc) irj = nsjac1(j3,jc) - DO jn=1,jpmsp - ij(jn)=njcoth(irj,jn) - lj(jn)=ij(jn) /= 0 - END DO - DO jn=1,2 - IF (lj(jn)) THEN - p = nonzero_map(ij(jn),ntro3) - spfj(:,p) = spfj(:,p) - & - prk(:,irj)*deriv(:,jc,1) ! *prkrat(jc) - END IF - END DO - DO jn=3,jpmsp - IF (lj(jn)) THEN - p = nonzero_map(ij(jn),ntro3) - spfj(:,p) = spfj(:,p) + & - prk(:,irj)*deriv(:,jc,1) ! *prkrat(jc) + DO jn = 1, jpmsp + ij = njcoth(irj,jn) + IF (ij /= 0) THEN + p = nonzero_map(ij,ntro3) + IF (jn < 3) THEN + spfj(:,p) = spfj(:,p) - prk(:,irj)*deriv(:,jc,1) + ELSE + spfj(:,p) = spfj(:,p) + prk(:,irj)*deriv(:,jc,1) + END IF END IF END DO END DO @@ -589,25 +557,23 @@ SUBROUTINE spfuljac(n_points,cdt,min_pivot,nonzero_map,spfj) ! 4. Add deposition terms to Jacobian diagonal. ! --- ---------- ----- -- -------- --------- ! -IF ( (ndepw /= 0) .OR. (ndepd /= 0) ) THEN +IF ((ndepw /= 0) .OR. (ndepd /= 0)) THEN DO js = 1, jpspec ifamd = moffam(js) itrd = madvtr(js) ityped = ctype(js) ! - IF ( ifamd /= 0 ) THEN - DO jl=1,n_points - IF ((ityped == jpfm) .OR. ((ityped == jpif) .AND. & - linfam(jl,itrd))) THEN - p = nonzero_map(ifamd,ifamd) - spfj(jl,p)=spfj(jl,p) & - - nodd(js)*(dpd(jl,js)+dpw(jl,js))*y(jl,js) + IF (ifamd /= 0) THEN + p = nonzero_map(ifamd,ifamd) + DO jl = 1, n_points + IF ((ityped == jpfm) .OR. ((ityped == jpif) .AND. linfam(jl,itrd))) THEN + spfj(jl,p) = spfj(jl,p) - nodd(js)*(dpd(jl,js) + dpw(jl,js))*y(jl,js) END IF END DO END IF - IF ( itrd /= 0 ) THEN + IF (itrd /= 0) THEN p = nonzero_map(itrd,itrd) - spfj(:,p) = spfj(:,p) - (dpd(:,js)+dpw(:,js))*y(:,js) + spfj(:,p) = spfj(:,p) - (dpd(:,js) + dpw(:,js))*y(:,js) END IF END DO END IF @@ -616,8 +582,8 @@ SUBROUTINE spfuljac(n_points,cdt,min_pivot,nonzero_map,spfj) ! 5. Jacobian elements in final form ! -------- -------- -- ----- ---- ! -DO p=1,total - spfj(:,p)=spfj(:,p)/f(:,base_tracer(p)) ! filter f earlier! +DO p = 1, total + spfj(:,p) = spfj(:,p) / f(:,base_tracer(p)) ! filter f earlier! END DO ! @@ -627,10 +593,10 @@ END SUBROUTINE spfuljac !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -SUBROUTINE splinslv2(n_points,bb,xx,min_pivot,max_val, & - nonzero_map_unordered,modified_map,spfj) +SUBROUTINE splinslv2(n_points, bb, xx, min_pivot, max_val, & + nonzero_map_unordered, modified_map, spfj) -USE asad_mod, ONLY: jpcspf, reorder, spfjsize_max, total +USE asad_mod, ONLY: jpcspf, spfjsize_max, total USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook @@ -670,6 +636,7 @@ SUBROUTINE splinslv2(n_points,bb,xx,min_pivot,max_val, & ! Local variables INTEGER :: total1 ! total number of nonzero entries in Jacobian INTEGER :: kr +INTEGER :: jl INTEGER :: i INTEGER :: j INTEGER :: ikr @@ -693,16 +660,16 @@ SUBROUTINE splinslv2(n_points,bb,xx,min_pivot,max_val, & ! Filter sparse Jacobian #if defined(IBM_XL_FORTRAN) ! Version optimised for IBM by using the fsel IBM-only intrinsic -DO j=1,total - DO i=1,n_points - tmp=fsel(spfj(i,j)+max_val,spfj(i,j),-max_val) - spfj(i,j)=fsel(tmp-max_val,max_val,tmp) +DO j = 1, total + DO jl = 1, n_points + tmp = fsel(spfj(jl,j)+max_val, spfj(jl,j), -max_val) + spfj(jl,j) = fsel(tmp-max_val, max_val, tmp) END DO END DO #else -DO j=1,total - DO i=1,n_points - spfj(i,j)=MIN(MAX(spfj(i,j),-max_val),max_val) +DO j = 1, total + DO jl = 1, n_points + spfj(jl,j) = MIN(MAX(spfj(jl,j), -max_val), max_val) END DO END DO #endif @@ -712,20 +679,20 @@ SUBROUTINE splinslv2(n_points,bb,xx,min_pivot,max_val, & ! The L U factors are overwritten onto the original sparse Jacobian array total1 = total modified_map(:,:) = nonzero_map_unordered(:,:) -DO kr=1,jpcspf - pivot = spfj(:,modified_map(kr,kr)) +DO kr = 1, jpcspf + pivot(:) = spfj(:,modified_map(kr,kr)) WHERE (ABS(pivot) > min_pivot) - pivot = 1.0/pivot + pivot(:) = 1.0 / pivot ELSE WHERE - pivot = max_val + pivot(:) = max_val END WHERE ! PIVOT = 1./spfj(:,modified_map(kr,kr)) - DO i=kr+1,jpcspf + DO i = kr+1, jpcspf ikr = modified_map(i,kr) IF (ikr > 0) THEN kfact = spfj(:,ikr)*pivot spfj(:,ikr) = kfact - DO j=kr+1,jpcspf + DO j = kr+1, jpcspf krj = modified_map(kr,j) IF (krj > 0) THEN ij = modified_map(i,j) @@ -747,57 +714,23 @@ SUBROUTINE splinslv2(n_points,bb,xx,min_pivot,max_val, & ! Filter sparse Jacobian #if defined(IBM_XL_FORTRAN) ! Version optimised for IBM by using the fsel IBM-only intrinsic -DO j=1,total1 - DO i=1,n_points - tmp=fsel(spfj(i,j)+max_val,spfj(i,j),-max_val) - spfj(i,j)=fsel(tmp-max_val,max_val,tmp) +DO j = 1, total1 + DO jl = 1, n_points + tmp = fsel(spfj(jl,j)+max_val, spfj(jl,j), -max_val) + spfj(jl,j) = fsel(tmp-max_val, max_val, tmp) END DO END DO #else -DO j=1,total1 - DO i=1,n_points - spfj(i,j)=MIN(MAX(spfj(i,j),-max_val),max_val) +DO j = 1, total1 + DO jl = 1, n_points + spfj(jl,j) = MIN(MAX(spfj(jl,j), -max_val), max_val) END DO END DO #endif ! Section 2: Solve P A P' z = P b with P'z = x using L U z = P b -! -! manual inlining of lu_solve() -! -! Part (a) calculate right-hand side P b -DO i=1,jpcspf - bb1(:,i) = bb(:,reorder(i)) -END DO - -! Part (b) forward-substitution: find w where L w = P b -DO kr=1,jpcspf-1 - DO i=kr+1,jpcspf - ikr = modified_map(i,kr) - IF (ikr > 0) & - bb1(:,i) = bb1(:,i) - spfj(:,ikr) * bb1(:,kr) - END DO -END DO - -! Part (c) back-substitution: find z where U z = w -DO kr=jpcspf,1,-1 - xx1(:,kr) = bb1(:,kr) - DO j=kr+1,jpcspf - krj = modified_map(kr,j) - IF (krj > 0) & - xx1(:,kr) = xx1(:,kr) - spfj(:,krj) * xx1(:,j) - END DO - pivot = spfj(:,modified_map(kr,kr)) - WHERE (ABS(pivot) < min_pivot) pivot = min_pivot - xx1(:,kr) = MIN(MAX(xx1(:,kr),-max_val),max_val)/pivot -END DO - -! Part (d) determine x, apply transpose(P) to z, P'z = x -xx = 0.0 -DO j=1,jpcspf - xx(:,reorder(j)) = xx1(:,j) -END DO +CALL spresolv2(n_points, bb, xx, min_pivot, modified_map, spfj, max_val) IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN @@ -805,7 +738,7 @@ END SUBROUTINE splinslv2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -SUBROUTINE spresolv2(n_points,bb,xx,min_pivot,modified_map,spfj) +SUBROUTINE spresolv2(n_points, bb, xx, min_pivot, modified_map, spfj, max_val) ! This subroutine determines x where L U z = P b with P' z = x ! The L U factors are supplied to this routine and contained with spfj array. @@ -830,6 +763,9 @@ SUBROUTINE spresolv2(n_points,bb,xx,min_pivot,modified_map,spfj) INTEGER, INTENT(IN) :: modified_map(jpcspf,jpcspf) REAL, INTENT(IN) :: spfj(n_points,spfjsize_max) +! Maximum tolerated value for use in filtering step. Unused if negative +REAL, INTENT(IN) :: max_val + INTEGER :: kr INTEGER :: i @@ -850,32 +786,37 @@ SUBROUTINE spresolv2(n_points,bb,xx,min_pivot,modified_map,spfj) IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) ! Part (a) calculate right-hand side P b -DO i=1,jpcspf +DO i = 1, jpcspf bb1(:,i) = bb(:,reorder(i)) END DO ! Part (b) forward-substitution: find w where L w = P b -DO kr=1,jpcspf-1 - DO i=kr+1,jpcspf +DO kr = 1, jpcspf-1 + DO i = kr+1, jpcspf ikr = modified_map(i,kr) IF (ikr > 0) bb1(:,i) = bb1(:,i) - spfj(:,ikr)*bb1(:,kr) END DO END DO ! Part (c) back-substitution: find z where U z = w -DO kr=jpcspf,1,-1 +DO kr = jpcspf, 1, -1 xx1(:,kr) = bb1(:,kr) - DO j=kr+1,jpcspf + DO j = kr+1, jpcspf krj = modified_map(kr,j) IF (krj > 0) xx1(:,kr) = xx1(:,kr) - spfj(:,krj)*xx1(:,j) END DO - pivot = spfj(:,modified_map(kr,kr)) + pivot(:) = spfj(:,modified_map(kr,kr)) WHERE (ABS(pivot) < min_pivot) pivot = min_pivot - xx1(:,kr) = xx1(:,kr) / pivot + IF (max_val > 0) THEN + xx1(:,kr) = MIN(MAX(xx1(:,kr), -max_val), max_val) / pivot + ELSE + xx1(:,kr) = xx1(:,kr) / pivot + END IF END DO -xx = 0.0 -DO i=1,jpcspf +! Part (d) determine x, apply transpose(P) to z, P'z = x +xx(:,:) = 0.0 +DO i = 1, jpcspf xx(:,reorder(i)) = xx1(:,i) END DO diff --git a/src/science/core/chemistry/asad/asad_spimpmjp.F90 b/src/science/core/chemistry/asad/asad_spimpmjp.F90 index 22d64c6..b162acd 100644 --- a/src/science/core/chemistry/asad/asad_spimpmjp.F90 +++ b/src/science/core/chemistry/asad/asad_spimpmjp.F90 @@ -294,6 +294,8 @@ SUBROUTINE asad_spimpmjp(exit_code, ix, jy, nlev, n_points, location, & ! The maximum concentration allowed was previously f_max = 1.0/f_min REAL, PARAMETER :: f_max = 1.0e30 REAL :: f_min +! Negative threshold to turn off filtering in quasi-Newton mode +REAL, PARAMETER :: max_val = -1.0 REAL :: RelTol_residual_error REAL :: RelTol_error REAL :: rafmin @@ -557,7 +559,8 @@ SUBROUTINE asad_spimpmjp(exit_code, ix, jy, nlev, n_points, location, & /DOT_PRODUCT(delta_G(jl,:),delta_G(jl,:)) G_ftmp(jl,:) = G_f(jl,:)*(1.0 - coeff) END DO - CALL spresolv2(n_points,G_ftmp,f_incr,f_min,modified_map,spfj) + + CALL spresolv2(n_points,G_ftmp,f_incr,f_min,modified_map,spfj,max_val) f = f + f_incr ! remove negative values. Does not need to be done in diff --git a/src/science/core/chemistry/asad/asad_steady.F90 b/src/science/core/chemistry/asad/asad_steady.F90 index 139ec90..a4a1e10 100644 --- a/src/science/core/chemistry/asad/asad_steady.F90 +++ b/src/science/core/chemistry/asad/asad_steady.F90 @@ -93,8 +93,6 @@ SUBROUTINE asad_steady( kl ) INTEGER, PARAMETER :: n_ho2=3 INTEGER, PARAMETER :: n_no=4 -INTEGER :: jl1 -INTEGER :: jl2 INTEGER :: jr INTEGER :: ix INTEGER :: i @@ -117,14 +115,6 @@ SUBROUTINE asad_steady( kl ) ! ! Set up loops correctly IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) -jl1 = 1 -jl2 = kl -IF (kl == 1) THEN - ! Previously the line below was jl1 = jl1+jlst, but jlst was never - ! set and so it has been removed. - jl1 = jl1 - jl2 = jl1 -END IF ! Initialise DERIV this timestep, first done in ASAD_INIT deriv(:,:,:) = 1.0 @@ -141,222 +131,190 @@ SUBROUTINE asad_steady( kl ) DO jr = 1,nsspt(ix) i = nsspi(ix,jr) IF (i <= nuni) THEN - ssnum(jl1:jl2) = ssnum(jl1:jl2) + & - rk(jl1:jl2,i)*y(jl1:jl2,nspi(i,1)) + ssnum(1:kl) = ssnum(1:kl) + rk(1:kl,i)*y(1:kl,nspi(i,1)) IF ((ix < 5) .AND. (nspi(i,1) == nspo3 )) & - ! add terms to derivative for d(j[O3])/d[O3] = j_o3 - dssnum(jl1:jl2,n_o3) = dssnum(jl1:jl2,n_o3) + & - rk(jl1:jl2,i) + ! add terms to derivative for d(j[O3])/d[O3] = j_o3 + dssnum(1:kl,n_o3) = dssnum(1:kl,n_o3) + rk(1:kl,i) + IF ((ix < 5) .AND. (nspi(i,1) == nspno )) & - ! add terms to derivative for d(j[NO])/d[NO] = j_no - dssnum(jl1:jl2,n_no) = dssnum(jl1:jl2,n_no) + & - rk(jl1:jl2,i) + ! add terms to derivative for d(j[NO])/d[NO] = j_no + dssnum(1:kl,n_no) = dssnum(1:kl,n_no) + rk(1:kl,i) + ELSE - ssnum(jl1:jl2) = ssnum(jl1:jl2) + & - rk(jl1:jl2,i)*y(jl1:jl2,nspi(i,1))*y(jl1:jl2,nspi(i,2)) + ssnum(1:kl) = ssnum(1:kl) + rk(1:kl,i)*y(1:kl,nspi(i,1))*y(1:kl,nspi(i,2)) IF (ix < 5) THEN ! add terms for derivative w.r.t. ozone. IF (nspi(i,1) == nspo1d) & - dssnum(jl1:jl2,n_o3) = dssnum(jl1:jl2,n_o3) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,2))*deriv(jl1:jl2,nss_o1d,n_o3) + dssnum(1:kl,n_o3) = dssnum(1:kl,n_o3) + & + rk(1:kl,i)*y(1:kl,nspi(i,2))*deriv(1:kl,nss_o1d,n_o3) IF (nspi(i,2) == nspo1d) & - dssnum(jl1:jl2,n_o3) = dssnum(jl1:jl2,n_o3) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,1))*deriv(jl1:jl2,nss_o1d,n_o3) + dssnum(1:kl,n_o3) = dssnum(1:kl,n_o3) + & + rk(1:kl,i)*y(1:kl,nspi(i,1))*deriv(1:kl,nss_o1d,n_o3) IF (nspi(i,1) == nspo3p) & - dssnum(jl1:jl2,n_o3) = dssnum(jl1:jl2,n_o3) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,2))*deriv(jl1:jl2,nss_o3p,n_o3) + dssnum(1:kl,n_o3) = dssnum(1:kl,n_o3) + & + rk(1:kl,i)*y(1:kl,nspi(i,2))*deriv(1:kl,nss_o3p,n_o3) IF (nspi(i,2) == nspo3p) & - dssnum(jl1:jl2,n_o3) = dssnum(jl1:jl2,n_o3) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,1))*deriv(jl1:jl2,nss_o3p,n_o3) + dssnum(1:kl,n_o3) = dssnum(1:kl,n_o3) + & + rk(1:kl,i)*y(1:kl,nspi(i,1))*deriv(1:kl,nss_o3p,n_o3) IF (nspi(i,1) == nspo3) & - dssnum(jl1:jl2,n_o3) = dssnum(jl1:jl2,n_o3) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,2)) + dssnum(1:kl,n_o3) = dssnum(1:kl,n_o3) + rk(1:kl,i)*y(1:kl,nspi(i,2)) IF (nspi(i,2) == nspo3) & - dssnum(jl1:jl2,n_o3) = dssnum(jl1:jl2,n_o3) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,1)) + dssnum(1:kl,n_o3) = dssnum(1:kl,n_o3) + rk(1:kl,i)*y(1:kl,nspi(i,1)) ! add terms for derivative w.r.t OH IF (nspi(i,1) == nspo3p) & - ! add terms to derivates for d(a[A][B]) - dssnum(jl1:jl2,n_oh) = dssnum(jl1:jl2,n_oh) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,2))*deriv(jl1:jl2,nss_o3p,n_oh) + ! add terms to derivates for d(a[A][B]) + dssnum(1:kl,n_oh) = dssnum(1:kl,n_oh) + & + rk(1:kl,i)*y(1:kl,nspi(i,2))*deriv(1:kl,nss_o3p,n_oh) IF (nspi(i,2) == nspo3p) & - ! add terms to derivates for d(a[O2][O1D])/d[O3] and b[N2][O1D] - dssnum(jl1:jl2,n_oh) = dssnum(jl1:jl2,n_oh) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,1))*deriv(jl1:jl2,nss_o3p,n_oh) + ! add terms to derivates for d(a[O2][O1D])/d[O3] and b[N2][O1D] + dssnum(1:kl,n_oh) = dssnum(1:kl,n_oh) + & + rk(1:kl,i)*y(1:kl,nspi(i,1))*deriv(1:kl,nss_o3p,n_oh) IF (nspi(i,1) == nspoh) & - ! add terms to derivates for d(a[O2][O1D])/d[O3] and b[N2][O1D] - dssnum(jl1:jl2,n_oh) = dssnum(jl1:jl2,n_oh) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,2)) + ! add terms to derivates for d(a[O2][O1D])/d[O3] and b[N2][O1D] + dssnum(1:kl,n_oh) = dssnum(1:kl,n_oh) + rk(1:kl,i)*y(1:kl,nspi(i,2)) IF (nspi(i,2) == nspoh) & - ! add terms to derivates for d(a[O2][O1D])/d[O3] and b[N2][O1D] - dssnum(jl1:jl2,n_oh) = dssnum(jl1:jl2,n_oh) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,1)) + ! add terms to derivates for d(a[O2][O1D])/d[O3] and b[N2][O1D] + dssnum(1:kl,n_oh) = dssnum(1:kl,n_oh) + rk(1:kl,i)*y(1:kl,nspi(i,1)) ! add terms for derivative w.r.t HO2 IF (nspi(i,1) == nspo3p) & - dssnum(jl1:jl2,n_ho2) = dssnum(jl1:jl2,n_ho2) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,2))*deriv(jl1:jl2,nss_o3p,n_ho2) + dssnum(1:kl,n_ho2) = dssnum(1:kl,n_ho2) + & + rk(1:kl,i)*y(1:kl,nspi(i,2))*deriv(1:kl,nss_o3p,n_ho2) IF (nspi(i,2) == nspo3p) & - dssnum(jl1:jl2,n_ho2) = dssnum(jl1:jl2,n_ho2) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,1))*deriv(jl1:jl2,nss_o3p,n_ho2) + dssnum(1:kl,n_ho2) = dssnum(1:kl,n_ho2) + & + rk(1:kl,i)*y(1:kl,nspi(i,1))*deriv(1:kl,nss_o3p,n_ho2) IF (nspi(i,1) == nspho2) & - dssnum(jl1:jl2,n_ho2) = dssnum(jl1:jl2,n_ho2) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,2)) + dssnum(1:kl,n_ho2) = dssnum(1:kl,n_ho2) + rk(1:kl,i)*y(1:kl,nspi(i,2)) IF (nspi(i,2) == nspho2) & - dssnum(jl1:jl2,n_ho2) = dssnum(jl1:jl2,n_ho2) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,1)) + dssnum(1:kl,n_ho2) = dssnum(1:kl,n_ho2) + rk(1:kl,i)*y(1:kl,nspi(i,1)) ! add terms for derivative w.r.t NO IF (nspi(i,1) == nspno) & - dssnum(jl1:jl2,n_no) = dssnum(jl1:jl2,n_no) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,2)) + dssnum(1:kl,n_no) = dssnum(1:kl,n_no) + rk(1:kl,i)*y(1:kl,nspi(i,2)) IF (nspi(i,2) == nspno) & - dssnum(jl1:jl2,n_no) = dssnum(jl1:jl2,n_no) + & - rk(jl1:jl2,i) & - *y(jl1:jl2,nspi(i,1)) + dssnum(1:kl,n_no) = dssnum(1:kl,n_no) + rk(1:kl,i)*y(1:kl,nspi(i,1)) END IF END IF - END DO ! jr + END DO ! jr ! ! Destruction terms DO jr = 1,nssrt(ix) i = nssri(ix,jr) j = nssrx(ix,jr) IF (i <= nuni) THEN - ssden(jl1:jl2) = ssden(jl1:jl2) + rk(jl1:jl2,i) + ssden(1:kl) = ssden(1:kl) + rk(1:kl,i) ELSE - ssden(jl1:jl2) = ssden(jl1:jl2) + & - rk(jl1:jl2,i) * y(jl1:jl2,nspi(i,j)) + ssden(1:kl) = ssden(1:kl) + rk(1:kl,i) * y(1:kl,nspi(i,j)) IF (ix < 5) THEN IF (nspi(i,j) == nspo3 ) & - dssden(jl1:jl2,n_o3 ) = dssden(jl1:jl2,n_o3 ) + & - rk(jl1:jl2,i) + dssden(1:kl,n_o3 ) = dssden(1:kl,n_o3 ) + rk(1:kl,i) IF (nspi(i,j) == nspoh ) & - dssden(jl1:jl2,n_oh ) = dssden(jl1:jl2,n_oh ) + & - rk(jl1:jl2,i) + dssden(1:kl,n_oh ) = dssden(1:kl,n_oh ) + rk(1:kl,i) IF (nspi(i,j) == nspho2) & - dssden(jl1:jl2,n_ho2) = dssden(jl1:jl2,n_ho2) + & - rk(jl1:jl2,i) + dssden(1:kl,n_ho2) = dssden(1:kl,n_ho2) + rk(1:kl,i) IF (nspi(i,j) == nspno ) & - dssden(jl1:jl2,n_no ) = dssden(jl1:jl2,n_no ) + & - rk(jl1:jl2,i) + dssden(1:kl,n_no ) = dssden(1:kl,n_no ) + rk(1:kl,i) END IF END IF - END DO ! jr + END DO ! jr ! ! Steady state and derivatives of steady state - y(jl1:jl2,nssi(ix)) = ssnum(jl1:jl2)/ssden(jl1:jl2) + y(1:kl,nssi(ix)) = ssnum(1:kl)/ssden(1:kl) IF (ix < 5) THEN DO jr =1,4 - deriv(jl1:jl2,ix,jr) = & - (ssden(jl1:jl2)*dssnum(jl1:jl2,jr) - & - ssnum(jl1:jl2)*dssden(jl1:jl2,jr))/ & - (ssden(jl1:jl2) * ssden(jl1:jl2)) - END DO ! jr + deriv(1:kl,ix,jr) = & + (ssden(1:kl)*dssnum(1:kl,jr) - ssnum(1:kl)*dssden(1:kl,jr)) / & + (ssden(1:kl)*ssden(1:kl)) + END DO ! jr END IF -END DO ! ix +END DO ! ix ! rescale deriv to mean [O3]/[O] * d[O]/d[O3], where [O] = [O(1D)] or [O(3P)] ! for O(1D), and O(3P), N and H when these are SS species -WHERE (y(jl1:jl2,nspo1d) > peps) - deriv(jl1:jl2,nss_o1d,n_o3) = deriv(jl1:jl2,nss_o1d,n_o3 )* & - y(jl1:jl2,nspo3 )/y(jl1:jl2,nspo1d) - deriv(jl1:jl2,nss_o1d,n_oh) = deriv(jl1:jl2,nss_o1d,n_oh )* & - y(jl1:jl2,nspoh )/y(jl1:jl2,nspo1d) - deriv(jl1:jl2,nss_o1d,n_ho2)= deriv(jl1:jl2,nss_o1d,n_ho2)* & - y(jl1:jl2,nspho2)/y(jl1:jl2,nspo1d) - deriv(jl1:jl2,nss_o1d,n_no )= deriv(jl1:jl2,nss_o1d,n_no )* & - y(jl1:jl2,nspno )/y(jl1:jl2,nspo1d) +WHERE (y(1:kl,nspo1d) > peps) + deriv(1:kl,nss_o1d,n_o3) = deriv(1:kl,nss_o1d,n_o3 )* & + y(1:kl,nspo3 )/y(1:kl,nspo1d) + deriv(1:kl,nss_o1d,n_oh) = deriv(1:kl,nss_o1d,n_oh )* & + y(1:kl,nspoh )/y(1:kl,nspo1d) + deriv(1:kl,nss_o1d,n_ho2)= deriv(1:kl,nss_o1d,n_ho2)* & + y(1:kl,nspho2)/y(1:kl,nspo1d) + deriv(1:kl,nss_o1d,n_no )= deriv(1:kl,nss_o1d,n_no )* & + y(1:kl,nspno )/y(1:kl,nspo1d) ELSE WHERE - deriv(jl1:jl2,nss_o1d,n_o3) = 1.0 - deriv(jl1:jl2,nss_o1d,n_oh) = 1.0 - deriv(jl1:jl2,nss_o1d,n_ho2) = 1.0 - deriv(jl1:jl2,nss_o1d,n_no) = 1.0 + deriv(1:kl,nss_o1d,n_o3) = 1.0 + deriv(1:kl,nss_o1d,n_oh) = 1.0 + deriv(1:kl,nss_o1d,n_ho2) = 1.0 + deriv(1:kl,nss_o1d,n_no) = 1.0 END WHERE IF (o3p_in_ss) THEN - WHERE (y(jl1:jl2,nspo3p) > peps) - deriv(jl1:jl2,nss_o3p,n_o3 )= deriv(jl1:jl2,nss_o3p,n_o3 )* & - y(jl1:jl2,nspo3 )/y(jl1:jl2,nspo3p) - deriv(jl1:jl2,nss_o3p,n_oh )= deriv(jl1:jl2,nss_o3p,n_oh )* & - y(jl1:jl2,nspoh )/y(jl1:jl2,nspo3p) - deriv(jl1:jl2,nss_o3p,n_ho2)= deriv(jl1:jl2,nss_o3p,n_ho2)* & - y(jl1:jl2,nspho2)/y(jl1:jl2,nspo3p) - deriv(jl1:jl2,nss_o3p,n_no )= deriv(jl1:jl2,nss_o3p,n_no )* & - y(jl1:jl2,nspno )/y(jl1:jl2,nspo3p) + WHERE (y(1:kl,nspo3p) > peps) + deriv(1:kl,nss_o3p,n_o3 )= deriv(1:kl,nss_o3p,n_o3 )* & + y(1:kl,nspo3 )/y(1:kl,nspo3p) + deriv(1:kl,nss_o3p,n_oh )= deriv(1:kl,nss_o3p,n_oh )* & + y(1:kl,nspoh )/y(1:kl,nspo3p) + deriv(1:kl,nss_o3p,n_ho2)= deriv(1:kl,nss_o3p,n_ho2)* & + y(1:kl,nspho2)/y(1:kl,nspo3p) + deriv(1:kl,nss_o3p,n_no )= deriv(1:kl,nss_o3p,n_no )* & + y(1:kl,nspno )/y(1:kl,nspo3p) ELSE WHERE - deriv(jl1:jl2,nss_o3p,n_o3) = 1.0 - deriv(jl1:jl2,nss_o3p,n_oh) = 1.0 - deriv(jl1:jl2,nss_o3p,n_ho2) = 1.0 - deriv(jl1:jl2,nss_o3p,n_no) = 1.0 + deriv(1:kl,nss_o3p,n_o3) = 1.0 + deriv(1:kl,nss_o3p,n_oh) = 1.0 + deriv(1:kl,nss_o3p,n_ho2) = 1.0 + deriv(1:kl,nss_o3p,n_no) = 1.0 END WHERE END IF IF (n_in_ss) THEN - WHERE (y(jl1:jl2,nspn ) > peps) - deriv(jl1:jl2,nss_n,n_o3) = deriv(jl1:jl2,nss_n,n_o3 )* & - y(jl1:jl2,nspo3 )/y(jl1:jl2,nspn) - deriv(jl1:jl2,nss_n,n_oh) = deriv(jl1:jl2,nss_n,n_oh )* & - y(jl1:jl2,nspoh )/y(jl1:jl2,nspn) - deriv(jl1:jl2,nss_n,n_ho2)= deriv(jl1:jl2,nss_n,n_ho2)* & - y(jl1:jl2,nspho2)/y(jl1:jl2,nspn) - deriv(jl1:jl2,nss_n,n_no )= deriv(jl1:jl2,nss_n,n_no )* & - y(jl1:jl2,nspno )/y(jl1:jl2,nspn) + WHERE (y(1:kl,nspn ) > peps) + deriv(1:kl,nss_n,n_o3) = deriv(1:kl,nss_n,n_o3 )* & + y(1:kl,nspo3 )/y(1:kl,nspn) + deriv(1:kl,nss_n,n_oh) = deriv(1:kl,nss_n,n_oh )* & + y(1:kl,nspoh )/y(1:kl,nspn) + deriv(1:kl,nss_n,n_ho2)= deriv(1:kl,nss_n,n_ho2)* & + y(1:kl,nspho2)/y(1:kl,nspn) + deriv(1:kl,nss_n,n_no )= deriv(1:kl,nss_n,n_no )* & + y(1:kl,nspno )/y(1:kl,nspn) ELSE WHERE - deriv(jl1:jl2,nss_n,n_o3) = 1.0 - deriv(jl1:jl2,nss_n,n_oh) = 1.0 - deriv(jl1:jl2,nss_n,n_ho2) = 1.0 - deriv(jl1:jl2,nss_n,n_no) = 1.0 + deriv(1:kl,nss_n,n_o3) = 1.0 + deriv(1:kl,nss_n,n_oh) = 1.0 + deriv(1:kl,nss_n,n_ho2) = 1.0 + deriv(1:kl,nss_n,n_no) = 1.0 END WHERE END IF IF (h_in_ss) THEN - WHERE (y(jl1:jl2,nsph ) > peps) - deriv(jl1:jl2,nss_h,n_o3) = deriv(jl1:jl2,nss_h,n_o3 )* & - y(jl1:jl2,nspo3 )/y(jl1:jl2,nsph ) - deriv(jl1:jl2,nss_h,n_oh) = deriv(jl1:jl2,nss_h,n_oh )* & - y(jl1:jl2,nspoh )/y(jl1:jl2,nsph ) - deriv(jl1:jl2,nss_h,n_ho2)= deriv(jl1:jl2,nss_h,n_ho2)* & - y(jl1:jl2,nspho2)/y(jl1:jl2,nsph ) - deriv(jl1:jl2,nss_h,n_no )= deriv(jl1:jl2,nss_h,n_no )* & - y(jl1:jl2,nspno )/y(jl1:jl2,nsph ) + WHERE (y(1:kl,nsph ) > peps) + deriv(1:kl,nss_h,n_o3) = deriv(1:kl,nss_h,n_o3 )* & + y(1:kl,nspo3 )/y(1:kl,nsph ) + deriv(1:kl,nss_h,n_oh) = deriv(1:kl,nss_h,n_oh )* & + y(1:kl,nspoh )/y(1:kl,nsph ) + deriv(1:kl,nss_h,n_ho2)= deriv(1:kl,nss_h,n_ho2)* & + y(1:kl,nspho2)/y(1:kl,nsph ) + deriv(1:kl,nss_h,n_no )= deriv(1:kl,nss_h,n_no )* & + y(1:kl,nspno )/y(1:kl,nsph ) ELSE WHERE - deriv(jl1:jl2,nss_h,n_o3) = 1.0 - deriv(jl1:jl2,nss_h,n_oh) = 1.0 - deriv(jl1:jl2,nss_h,n_ho2) = 1.0 - deriv(jl1:jl2,nss_h,n_no) = 1.0 + deriv(1:kl,nss_h,n_o3) = 1.0 + deriv(1:kl,nss_h,n_oh) = 1.0 + deriv(1:kl,nss_h,n_ho2) = 1.0 + deriv(1:kl,nss_h,n_no) = 1.0 END WHERE END IF diff --git a/src/science/core/chemistry/ukca_chemistry_ctl.F90 b/src/science/core/chemistry/ukca_chemistry_ctl.F90 index a8ba541..192e09a 100644 --- a/src/science/core/chemistry/ukca_chemistry_ctl.F90 +++ b/src/science/core/chemistry/ukca_chemistry_ctl.F90 @@ -196,7 +196,6 @@ SUBROUTINE ukca_chemistry_ctl( & REAL, ALLOCATABLE :: ystore(:) ! array for H2SO4 when updated in MODE REAL :: zftr(theta_field_size,jpcspf) ! 1-D array of chemically active species ! including RO2 species, in VMR -REAL :: cdot(theta_field_size,jpcspf) ! 1-D chem. tendency REAL :: zq(theta_field_size) ! 1-D water vapour vmr REAL :: co2_1d(theta_field_size) ! 1-D CO2 REAL :: zprt1d(theta_field_size,jppj) ! 1-D photolysis rates for ASAD @@ -239,8 +238,8 @@ SUBROUTINE ukca_chemistry_ctl( & END DO !$OMP PARALLEL DEFAULT(NONE) & -!$OMP PRIVATE(cdot, ddmask, errcode, ierr, jna, jro2, js, jspf, jtr, k, kcs, & -!$OMP kce, l, rc_het, ystore, zdryrt2, zftr, zprt1d, zq, co2_1d) & +!$OMP PRIVATE(ddmask, errcode, ierr, jna, jro2, js, jspf, jtr, k, kcs, kce, l, & +!$OMP rc_het, ystore, zdryrt2, zftr, zprt1d, zq, co2_1d) & !$OMP SHARED(advt, atm_cf2cl2_mol, atm_cfcl3_mol, atm_ch4_mol, & !$OMP atm_co_mol, atm_h2_mol, atm_mebr_mol, atm_n2o_mol, avogadro, & !$OMP c_species, c_na_species, cloud_frac, cmessage, & @@ -413,8 +412,7 @@ SUBROUTINE ukca_chemistry_ctl( & END IF END IF - CALL asad_cdrive(cdot, & - zftr, & + CALL asad_cdrive(zftr, & pres(kcs:kce), & temp(kcs:kce), & zq, & diff --git a/src/science/core/chemistry/ukca_chemistry_ctl_col_mod.F90 b/src/science/core/chemistry/ukca_chemistry_ctl_col_mod.F90 index f184e92..820d81f 100644 --- a/src/science/core/chemistry/ukca_chemistry_ctl_col_mod.F90 +++ b/src/science/core/chemistry/ukca_chemistry_ctl_col_mod.F90 @@ -207,7 +207,6 @@ SUBROUTINE ukca_chemistry_ctl_col( & REAL :: zt (model_levels) ! 1-D temperature REAL :: zclw(model_levels) ! 1-D cloud liquid water REAL :: zfcloud(model_levels) ! 1-D cloud fraction -REAL :: cdot(model_levels,jpcspf) ! 1-D chem. tendency REAL :: zq(model_levels) ! 1-D water vapour vmr REAL :: co2_1d(model_levels) ! 1-D CO2 vmr REAL :: zprt1d(model_levels,jppj) ! 1-D photolysis rates for ASAD @@ -280,8 +279,7 @@ SUBROUTINE ukca_chemistry_ctl_col( & ! Model levels loop !$OMP PARALLEL DEFAULT(NONE) & -!$OMP PRIVATE(cdot, cmessage, errcode, i, ierr, & -!$OMP j, js, l, rc_het, stratflag, & +!$OMP PRIVATE(cmessage, errcode, i, ierr, j, js, l, rc_het, stratflag, & !$OMP ystore, zclw, zdryrt2, zfcloud, zftr, have_nat1d, & !$OMP zp, zprt1d, zq, zt, co2_1d, zwetrt2, & !$OMP kcs, kce, chunk_size, dpd_full, dpw_full, & @@ -495,8 +493,7 @@ SUBROUTINE ukca_chemistry_ctl_col( & END IF ! Call asad_cdrive with segmented arrays - CALL asad_cdrive(cdot(kcs:kce,:), & - zftr(kcs:kce,:), & + CALL asad_cdrive(zftr(kcs:kce,:), & zp(kcs:kce), & zt(kcs:kce), & zq(kcs:kce), & diff --git a/src/science/core/chemistry/ukca_chemistry_ctl_full_mod.F90 b/src/science/core/chemistry/ukca_chemistry_ctl_full_mod.F90 index c42f586..e5c99bd 100644 --- a/src/science/core/chemistry/ukca_chemistry_ctl_full_mod.F90 +++ b/src/science/core/chemistry/ukca_chemistry_ctl_full_mod.F90 @@ -189,7 +189,6 @@ SUBROUTINE ukca_chemistry_ctl_full( & REAL, ALLOCATABLE :: ystore(:) ! array for H2SO4 when updated in MODE REAL :: zftr(tot_n_pnts,jpcspf) ! 1-D array of chemically active species ! including RO2 species, in VMR -REAL :: cdot(tot_n_pnts,jpcspf) ! 1-D chem. tendency REAL :: zq(tot_n_pnts) ! 1-D water vapour vmr REAL :: co2_1d(tot_n_pnts) ! 1-D CO2 REAL :: zprt1d(tot_n_pnts,jppj) ! 1-D photolysis rates for ASAD @@ -374,8 +373,8 @@ SUBROUTINE ukca_chemistry_ctl_full( & ! Note: unlike in the other chemistry_ctl variants, here we call asad_cdrive ! once, passing in all data points from the current MPI process. -CALL asad_cdrive(cdot, zftr, pres, temp, zq, co2_1d, cloud_frac, qcl, ix, jy, & - k, zdryrt2, zwetrt, rc_het, zprt1d, tot_n_pnts, have_nat, & +CALL asad_cdrive(zftr, pres, temp, zq, co2_1d, cloud_frac, qcl, ix, jy, k, & + zdryrt2, zwetrt, rc_het, zprt1d, tot_n_pnts, have_nat, & stratflag, H_plus) IF (ukca_config%l_ukca_het_psc) THEN From eda813c71eaf3d86b8e7dfcc5e562c9ad1e48fe2 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Wed, 17 Dec 2025 14:49:57 +0000 Subject: [PATCH 2/7] Migrate r7184_further_chem_refactor over from FCM --- src/science/core/chemistry/asad/asad_mod.F90 | 3 - .../core/chemistry/asad/asad_sparse_vars.F90 | 153 +++++++++--------- .../core/chemistry/asad/asad_spimpmjp.F90 | 12 +- .../core/chemistry/asad/asad_spmjpdriv.F90 | 1 + 4 files changed, 80 insertions(+), 89 deletions(-) diff --git a/src/science/core/chemistry/asad/asad_mod.F90 b/src/science/core/chemistry/asad/asad_mod.F90 index 60c8e1a..6af884c 100644 --- a/src/science/core/chemistry/asad/asad_mod.F90 +++ b/src/science/core/chemistry/asad/asad_mod.F90 @@ -166,7 +166,6 @@ MODULE asad_mod REAL, ALLOCATABLE :: frph(:) ! fractional product array (het) REAL, ALLOCATABLE :: frpx(:) ! fractional product array (total) ! sparse algebra -INTEGER, ALLOCATABLE :: nonzero_map_unordered(:,:) ! Map of nonzero entries INTEGER, ALLOCATABLE :: modified_map(:,:) ! modified map (after decomposition) INTEGER, ALLOCATABLE :: nonzero_map(:,:) ! Map of nonzero entries, before reordering @@ -542,8 +541,6 @@ SUBROUTINE asad_mod_init(n_points) njcoth(:,:) = 0 IF (method == int_method_NR) THEN - IF (.NOT. ALLOCATED(nonzero_map_unordered)) & - ALLOCATE(nonzero_map_unordered(jpcspf, jpcspf)) IF (.NOT. ALLOCATED(modified_map)) ALLOCATE(modified_map(jpcspf, jpcspf)) IF (.NOT. ALLOCATED(nonzero_map)) ALLOCATE(nonzero_map(jpcspf, jpcspf)) IF (.NOT. ALLOCATED(reorder)) ALLOCATE(reorder(jpcspf)) diff --git a/src/science/core/chemistry/asad/asad_sparse_vars.F90 b/src/science/core/chemistry/asad/asad_sparse_vars.F90 index 8bc69df..c1700d5 100644 --- a/src/science/core/chemistry/asad/asad_sparse_vars.F90 +++ b/src/science/core/chemistry/asad/asad_sparse_vars.F90 @@ -81,20 +81,21 @@ SUBROUTINE setup_spfuljac() ! ! 1. Create nonzero_map(1:jpcspf,1:jpcspf) and base_tracer(1:spfjsize_max) ! -! nonzero_map(1:jpcspf,1:jpcspf) is an integer indexing array such that -! spfj(nonzero_map(i,j)) = A(i,j) where A is dense Jacobian matrix and -! spfj is the compressed (sparse) storage of the dense Jacobian matrix. +! nonzero_map(1:jpcspf,1:jpcspf) is an integer indexing array such that, for +! a given species index jl, spfj(jl,nonzero_map(i,j)) = A(i,j) where A is +! the dense Jacobian matrix for species jl and spfj is the compressed +! (sparse) storage of the dense Jacobian matrix. ! -! 2. Create nonzero_map_unordered(1:jpcspf,1:jpcspf) +! 2. Create permuted_nonzero_map(1:jpcspf,1:jpcspf) ! -! nonzero_map_unordered(1:jpcspf,1:jpcspf) is an integer indexing array such -! that spfj(nonzero_map_unordered(i,j)) = (PAP')(i,j) where (PAP')(i,j) is +! permuted_nonzero_map(1:jpcspf,1:jpcspf) is an integer indexing array such +! that spfj(jl,permuted_nonzero_map(i,j)) = (PAP')(i,j) where (PAP')(i,j) is ! the (i,j) entry of the matrix matmul(P,matmul(A,transpose(P))) where P is ! the permutation matrix and A is the dense Jacobian. ! ! 3. Check number of nonzero entries in LU factorization of the dense Jacobian ! -! Gaussian-elimination is used to solve a linear equation in splinslv2 and +! Gaussian elimination is used to solve a linear equation in splinslv2 and ! the array spfj is used to hold the LU factorization. This section checks ! that when this LU factorization is done the spfj array is sufficiently ! large enough to hold the nonzero entries of the LU factorization. @@ -106,7 +107,7 @@ SUBROUTINE setup_spfuljac() USE asad_mod, ONLY: specf, frpx, jpcspf, jpfrpx, jpmsp, jpspec, & madvtr, modified_map, ndepd, ndepw, nfrpx, njcoth, nltrf, & - nmsjac, nmzjac, nonzero_map, nonzero_map_unordered, & + nmsjac, nmzjac, nonzero_map, & npdfr, nsjac1, nstst, ntabpd, ntrf, ntro3, nzjac1, & reorder, spfjsize_max, maxterms, maxfterms, & nposterms, nnegterms, nfracterms, posterms, negterms, & @@ -146,10 +147,9 @@ SUBROUTINE setup_spfuljac() INTEGER :: itemp1 INTEGER :: activity(jpcspf) -INTEGER, ALLOCATABLE :: permute(:,:) ! permutation matrix -INTEGER, ALLOCATABLE :: map(:,:) ! a matrix of 1's indicating where the - ! dense Jacobian A(i,j) has nonzero - ! entries +INTEGER :: permute(jpcspf,jpcspf) ! permutation matrix +INTEGER :: map(jpcspf,jpcspf) ! sparsity pattern of Jacobian +INTEGER :: permuted_nonzero_map(jpcspf,jpcspf) ! permuted map of nonzeros CHARACTER(LEN=errormessagelength) :: cmessage @@ -166,8 +166,6 @@ SUBROUTINE setup_spfuljac() ! base_tracer ! ------------------------------------------------------------------------------ -IF (.NOT. ALLOCATED(map)) ALLOCATE(map(jpcspf, jpcspf)) - map(:,:) = 0 DO i = 1, jpcspf map(i,i) = 1 @@ -183,10 +181,9 @@ SUBROUTINE setup_spfuljac() IF (ij /= 0) map(ij,itrcr) = 1 END DO ! - IF (npdfr(irj,1) /= 0) THEN - i1 = npdfr(irj,1) - i2 = npdfr(irj,2) - DO jn = i1, i2 + i1 = npdfr(irj,1) + IF (i1 /= 0) THEN + DO jn = i1, npdfr(irj,2) isx = ntabpd(jn,1) map(isx,itrcr) = 1 END DO @@ -241,15 +238,15 @@ SUBROUTINE setup_spfuljac() END DO base_tracer(:) = 0 -DO i = 1, jpcspf - DO j = 1, jpcspf - i1 = nonzero_map(j,i) - IF (i1 > 0) base_tracer(i1) = i +DO j = 1, jpcspf + DO i = 1, jpcspf + i1 = nonzero_map(i,j) + IF (i1 > 0) base_tracer(i1) = j END DO END DO ! ------------------------------------------------------------------------------ -! Section 2: Create nonzero_map_unordered (index array for PAP') +! Section 2: Create permuted_nonzero_map (index array for PAP') ! ------------------------------------------------------------------------------ ! Reorder species by their reactivity to minimize fill-in @@ -257,11 +254,9 @@ SUBROUTINE setup_spfuljac() activity(i) = SUM(map(:,i)) + SUM(map(i,:)) reorder(i) = i END DO - DO i = 1, jpcspf-1 DO j = i+1, jpcspf IF (activity(i) > activity(j)) THEN - ! exchange i and j tracers if i is more active than j. itemp1 = reorder(i) reorder(i) = reorder(j) reorder(j) = itemp1 @@ -280,22 +275,18 @@ SUBROUTINE setup_spfuljac() END IF END DO -! reorganize pointer variable to account for varying fill-in -IF (.NOT. ALLOCATED(permute)) ALLOCATE(permute(jpcspf,jpcspf)) - +! Determine permutation matrix P permute(:,:) = 0 DO i = 1, jpcspf permute(i,reorder(i)) = 1 END DO -! Calculate the index array nonzero_map_unordered such that -! spfj(nonzero_map_unordered(i,j)) = P*A*P'(i,j) -! where P is the permutation matrix and A is the dense Jacobian matrix -! with spfj the compressed storage (sparse) array for matrix A. -nonzero_map_unordered = MATMUL(MATMUL(permute, nonzero_map), TRANSPOSE(permute)) - -IF (ALLOCATED(permute)) DEALLOCATE(permute) -IF (ALLOCATED(map)) DEALLOCATE(map) +! Calculate the index array permuted_nonzero_map such that +! spfj(jl,permuted_nonzero_map(i,j)) = P*A*P'(i,j) +! for each chemical species jl where P is the permutation matrix and A is the +! dense Jacobian matrix with spfj the compressed storage (sparse) array for +! matrix A. +permuted_nonzero_map = MATMUL(MATMUL(permute, nonzero_map), TRANSPOSE(permute)) ! ------------------------------------------------------------------------------ ! Section 3: Check the number of nonzero elements in LU factorization @@ -304,7 +295,7 @@ SUBROUTINE setup_spfuljac() ! Calculate the number of nonzero matrix elements in the LU factorization of the ! array PAP' and check that it is less than spfjsize_max. total1 = total -modified_map(:,:) = nonzero_map_unordered(:,:) +modified_map(:,:) = permuted_nonzero_map(:,:) DO kr = 1, jpcspf DO i = kr+1, jpcspf ikr = modified_map(i,kr) @@ -313,8 +304,7 @@ SUBROUTINE setup_spfuljac() krj = modified_map(kr,j) IF (krj > 0) THEN ij = modified_map(i,j) - ! Distinguish whether matrix element is zero or not. If not, proceed - ! as in dense case. If it is, create new matrix element. + ! If the matrix entry is zero then create new matrix element. IF (ij <= 0) THEN total1 = total1 + 1 modified_map(i,j) = total1 @@ -325,7 +315,6 @@ SUBROUTINE setup_spfuljac() END DO END DO -! Perform error check outside of the loop to better suit GPU runs IF (total1 > spfjsize_max) THEN errcode = total1 WRITE(umMessage,'(A,2I4)') 'Total1 exceeded spfjsize_max: ', total1, & @@ -449,7 +438,7 @@ END SUBROUTINE setup_spfuljac !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -SUBROUTINE spfuljac(n_points, cdt, min_pivot, nonzero_map, spfj) +SUBROUTINE spfuljac(n_points, cdt, min_pivot, spfj) ! ! Routine to calculate the Jacobian in sparse format ! @@ -457,7 +446,7 @@ SUBROUTINE spfuljac(n_points, cdt, min_pivot, nonzero_map, spfj) jpspec, linfam, madvtr, moffam, ndepd, ndepw, njcoth, & nmsjac, nodd, nsjac1, nstst, ntro3, prk, spfjsize_max, & nposterms, nnegterms, nfracterms, posterms, negterms, & - fracterms, base_tracer, ffrac, y, total + fracterms, base_tracer, ffrac, y, total, nonzero_map USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook @@ -467,7 +456,6 @@ SUBROUTINE spfuljac(n_points, cdt, min_pivot, nonzero_map, spfj) INTEGER, INTENT(IN) :: n_points REAL, INTENT(IN) :: cdt REAL, INTENT(IN) :: min_pivot -INTEGER, INTENT(IN) :: nonzero_map(jpcspf,jpcspf) REAL, INTENT(OUT) :: spfj(n_points,spfjsize_max) ! Local variables @@ -561,9 +549,9 @@ SUBROUTINE spfuljac(n_points, cdt, min_pivot, nonzero_map, spfj) DO js = 1, jpspec ifamd = moffam(js) itrd = madvtr(js) - ityped = ctype(js) ! IF (ifamd /= 0) THEN + ityped = ctype(js) p = nonzero_map(ifamd,ifamd) DO jl = 1, n_points IF ((ityped == jpfm) .OR. ((ityped == jpif) .AND. linfam(jl,itrd))) THEN @@ -593,10 +581,9 @@ END SUBROUTINE spfuljac !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -SUBROUTINE splinslv2(n_points, bb, xx, min_pivot, max_val, & - nonzero_map_unordered, modified_map, spfj) +SUBROUTINE splinslv2(n_points, bb, xx, min_pivot, max_val, spfj) -USE asad_mod, ONLY: jpcspf, spfjsize_max, total +USE asad_mod, ONLY: jpcspf, modified_map, spfjsize_max, total USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook @@ -629,12 +616,9 @@ SUBROUTINE splinslv2(n_points, bb, xx, min_pivot, max_val, & REAL, INTENT(OUT) :: xx(n_points,jpcspf) REAL, INTENT(IN) :: min_pivot REAL, INTENT(IN) :: max_val -INTEGER, INTENT(IN) :: nonzero_map_unordered(jpcspf,jpcspf) -INTEGER, INTENT(IN OUT) :: modified_map(jpcspf,jpcspf) REAL, INTENT(IN OUT) :: spfj(1:n_points,1:spfjsize_max) ! Local variables -INTEGER :: total1 ! total number of nonzero entries in Jacobian INTEGER :: kr INTEGER :: jl INTEGER :: i @@ -643,11 +627,9 @@ SUBROUTINE splinslv2(n_points, bb, xx, min_pivot, max_val, & INTEGER :: krj INTEGER :: ij -REAL :: bb1(n_points,jpcspf) -REAL :: xx1(n_points,jpcspf) - +REAL :: multiplier(n_points) REAL :: pivot(n_points) -REAL :: kfact(n_points) +REAL :: lower(n_points) INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 @@ -657,6 +639,10 @@ SUBROUTINE splinslv2(n_points, bb, xx, min_pivot, max_val, & IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) +! NOTE: splinslv2 is always preceded by spfuljac, which computes spfj to be a +! sparse (compressed) representation of the Jacobian matrix for each +! species. + ! Filter sparse Jacobian #if defined(IBM_XL_FORTRAN) ! Version optimised for IBM by using the fsel IBM-only intrinsic @@ -676,35 +662,46 @@ SUBROUTINE splinslv2(n_points, bb, xx, min_pivot, max_val, & ! Section 1: Determine L U factors such that L U = P A P' -! The L U factors are overwritten onto the original sparse Jacobian array -total1 = total -modified_map(:,:) = nonzero_map_unordered(:,:) +! +! The L U factors are overwritten onto the original sparse Jacobian array. That +! is, for each species index jl, the upper triangular part of the matrix +! represented by spfj(jl,:) holds the U factor and the strictly lower triangular +! part holds the L factor below its unit diagonal. + +! Loop over rows for pivoting DO kr = 1, jpcspf - pivot(:) = spfj(:,modified_map(kr,kr)) - WHERE (ABS(pivot) > min_pivot) - pivot(:) = 1.0 / pivot + + ! The multipliers are the diagonal entries for each species + multiplier(:) = spfj(:,modified_map(kr,kr)) + + ! Determine the pivot for each species using a threshold + WHERE (ABS(multiplier) > min_pivot) + pivot(:) = 1.0 / multiplier ELSE WHERE pivot(:) = max_val END WHERE - ! PIVOT = 1./spfj(:,modified_map(kr,kr)) + + ! Loop over non-zero entries in the lower-diagonal part of column kr DO i = kr+1, jpcspf ikr = modified_map(i,kr) IF (ikr > 0) THEN - kfact = spfj(:,ikr)*pivot - spfj(:,ikr) = kfact + + ! Compute contribution to the lower diagonal part and store it in spfj + lower(:) = spfj(:,ikr)*pivot + spfj(:,ikr) = lower + + ! Loop over non-zero entries in the upper-diagonal part of row kr DO j = kr+1, jpcspf krj = modified_map(kr,j) IF (krj > 0) THEN + + ! Compute the contribution towards the upper triangular part and store + ! it in the upper-diagonal part of spfj + ! NOTE: We can rely on the entry being nonzero because the sparsity + ! pattern was pre-computed. ij = modified_map(i,j) - ! Distinguish whether matrix element is zero or not. If not, proceed - ! as in dense case. If it is, create new matrix element. - IF (ij > 0) THEN - spfj(:,ij) = spfj(:,ij) - kfact*spfj(:,krj) - ELSE - total1 = total1 + 1 - modified_map(i,j) = total1 - spfj(:,total1) = -kfact*spfj(:,krj) - END IF + spfj(:,ij) = spfj(:,ij) - lower*spfj(:,krj) + END IF END DO END IF @@ -714,14 +711,14 @@ SUBROUTINE splinslv2(n_points, bb, xx, min_pivot, max_val, & ! Filter sparse Jacobian #if defined(IBM_XL_FORTRAN) ! Version optimised for IBM by using the fsel IBM-only intrinsic -DO j = 1, total1 +DO j = 1, total DO jl = 1, n_points tmp = fsel(spfj(jl,j)+max_val, spfj(jl,j), -max_val) spfj(jl,j) = fsel(tmp-max_val, max_val, tmp) END DO END DO #else -DO j = 1, total1 +DO j = 1, total DO jl = 1, n_points spfj(jl,j) = MIN(MAX(spfj(jl,j), -max_val), max_val) END DO @@ -730,7 +727,7 @@ SUBROUTINE splinslv2(n_points, bb, xx, min_pivot, max_val, & ! Section 2: Solve P A P' z = P b with P'z = x using L U z = P b -CALL spresolv2(n_points, bb, xx, min_pivot, modified_map, spfj, max_val) +CALL spresolv2(n_points, bb, xx, min_pivot, spfj, max_val) IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN @@ -738,7 +735,7 @@ END SUBROUTINE splinslv2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -SUBROUTINE spresolv2(n_points, bb, xx, min_pivot, modified_map, spfj, max_val) +SUBROUTINE spresolv2(n_points, bb, xx, min_pivot, spfj, max_val) ! This subroutine determines x where L U z = P b with P' z = x ! The L U factors are supplied to this routine and contained with spfj array. @@ -749,7 +746,7 @@ SUBROUTINE spresolv2(n_points, bb, xx, min_pivot, modified_map, spfj, max_val) ! Part (c) back-substitution: find z where U z = w ! Part (d) determine x, apply transpose(P) to z, P'z = x -USE asad_mod, ONLY: jpcspf, reorder, spfjsize_max +USE asad_mod, ONLY: jpcspf, modified_map, reorder, spfjsize_max USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook @@ -760,13 +757,11 @@ SUBROUTINE spresolv2(n_points, bb, xx, min_pivot, modified_map, spfj, max_val) REAL, INTENT(IN) :: bb(n_points,jpcspf) REAL, INTENT(OUT) :: xx(n_points,jpcspf) REAL, INTENT(IN) :: min_pivot -INTEGER, INTENT(IN) :: modified_map(jpcspf,jpcspf) REAL, INTENT(IN) :: spfj(n_points,spfjsize_max) ! Maximum tolerated value for use in filtering step. Unused if negative REAL, INTENT(IN) :: max_val - INTEGER :: kr INTEGER :: i INTEGER :: j diff --git a/src/science/core/chemistry/asad/asad_spimpmjp.F90 b/src/science/core/chemistry/asad/asad_spimpmjp.F90 index b162acd..ded7860 100644 --- a/src/science/core/chemistry/asad/asad_spimpmjp.F90 +++ b/src/science/core/chemistry/asad/asad_spimpmjp.F90 @@ -255,8 +255,7 @@ SUBROUTINE asad_spimpmjp(exit_code, ix, jy, nlev, n_points, location, & solver_iter) USE asad_mod, ONLY: ptol, peps, cdt, f, fdot, nitnr, nstst, y, & - fj, nonzero_map, ltrig, jpcspf, spfj, & - modified_map, nonzero_map_unordered + fj, nonzero_map, ltrig, jpcspf, spfj USE asad_sparse_vars, ONLY: setup_spfuljac, spfuljac, spresolv2, splinslv2 USE ukca_config_specification_mod, ONLY: ukca_config USE yomhook, ONLY: lhook, dr_hook @@ -375,7 +374,7 @@ SUBROUTINE asad_spimpmjp(exit_code, ix, jy, nlev, n_points, location, & END IF !$OMP END CRITICAL (setup_jacobian_init) -CALL spfuljac(n_points,cdt,f_min,nonzero_map,spfj) +CALL spfuljac(n_points,cdt,f_min,spfj) ! Call forward Euler to make first guess, f_0 for f(t=n+1) (next timestep) CALL forward_euler(n_points, f, f_initial, f_min, nonzero_map, spfj) @@ -443,7 +442,7 @@ SUBROUTINE asad_spimpmjp(exit_code, ix, jy, nlev, n_points, location, & END DO END DO - CALL spfuljac(n_points,cdt,f_min,nonzero_map,spfj) + CALL spfuljac(n_points,cdt,f_min,spfj) IF (ltrig .AND. printstatus == PrStatus_Diag) THEN WRITE(umMessage,"('Iteration ',i4)") iter @@ -465,8 +464,7 @@ SUBROUTINE asad_spimpmjp(exit_code, ix, jy, nlev, n_points, location, & END DO END IF - CALL splinslv2(n_points,G_f,f_incr,f_min,f_max,nonzero_map_unordered, & - modified_map,spfj) + CALL splinslv2(n_points,G_f,f_incr,f_min,f_max,spfj) IF (ltrig .AND. printstatus == PrStatus_Diag) THEN DO jl=1,n_points @@ -560,7 +558,7 @@ SUBROUTINE asad_spimpmjp(exit_code, ix, jy, nlev, n_points, location, & G_ftmp(jl,:) = G_f(jl,:)*(1.0 - coeff) END DO - CALL spresolv2(n_points,G_ftmp,f_incr,f_min,modified_map,spfj,max_val) + CALL spresolv2(n_points,G_ftmp,f_incr,f_min,spfj,max_val) f = f + f_incr ! remove negative values. Does not need to be done in diff --git a/src/science/core/chemistry/asad/asad_spmjpdriv.F90 b/src/science/core/chemistry/asad/asad_spmjpdriv.F90 index 599b680..e380c23 100644 --- a/src/science/core/chemistry/asad/asad_spmjpdriv.F90 +++ b/src/science/core/chemistry/asad/asad_spmjpdriv.F90 @@ -117,6 +117,7 @@ SUBROUTINE asad_spmjpdriv(ix,jy,nlev,n_points,num_iter) USE errormessagelength_mod, ONLY: errormessagelength +! TODO: This import will become redundant USE asad_diffun_mod, ONLY: asad_diffun USE asad_spimpmjp_mod, ONLY: asad_spimpmjp USE asad_ftoy_mod, ONLY: asad_ftoy From 7c99b0d9f226c06f1a356fa3de0e6cf398f2aa4f Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Fri, 23 Jan 2026 10:55:58 +0000 Subject: [PATCH 3/7] Post-merge fix --- src/science/core/chemistry/asad/asad_mod.F90 | 3 +++ src/science/core/chemistry/asad/asad_sparse_vars.F90 | 3 +-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/science/core/chemistry/asad/asad_mod.F90 b/src/science/core/chemistry/asad/asad_mod.F90 index 0c6157d..5637920 100644 --- a/src/science/core/chemistry/asad/asad_mod.F90 +++ b/src/science/core/chemistry/asad/asad_mod.F90 @@ -166,6 +166,7 @@ MODULE asad_mod REAL, ALLOCATABLE :: frph(:) ! fractional product array (het) REAL, ALLOCATABLE :: frpx(:) ! fractional product array (total) ! sparse algebra +INTEGER, ALLOCATABLE :: permuted_nonzero_map(:,:) ! permuted map of nonzeros INTEGER, ALLOCATABLE :: modified_map(:,:) ! modified map (after decomposition) INTEGER, ALLOCATABLE :: nonzero_map(:,:) ! Map of nonzero entries, before reordering @@ -543,6 +544,8 @@ SUBROUTINE asad_mod_init(n_points) njcoth(:,:) = 0 IF (method == int_method_NR) THEN + IF (.NOT. ALLOCATED(permuted_nonzero_map)) & + ALLOCATE(permuted_nonzero_map(jpcspf, jpcspf)) IF (.NOT. ALLOCATED(modified_map)) ALLOCATE(modified_map(jpcspf, jpcspf)) IF (.NOT. ALLOCATED(nonzero_map)) ALLOCATE(nonzero_map(jpcspf, jpcspf)) IF (.NOT. ALLOCATED(reorder)) ALLOCATE(reorder(jpcspf)) diff --git a/src/science/core/chemistry/asad/asad_sparse_vars.F90 b/src/science/core/chemistry/asad/asad_sparse_vars.F90 index aa86431..5515e39 100644 --- a/src/science/core/chemistry/asad/asad_sparse_vars.F90 +++ b/src/science/core/chemistry/asad/asad_sparse_vars.F90 @@ -107,7 +107,7 @@ SUBROUTINE setup_spfuljac() USE asad_mod, ONLY: specf, frpx, jpcspf, jpfrpx, jpmsp, jpspec, & madvtr, modified_map, ndepd, ndepw, nfrpx, njcoth, nltrf, & - nmsjac, nmzjac, nonzero_map, & + nmsjac, nmzjac, nonzero_map, permuted_nonzero_map, & npdfr, nsjac1, nstst, ntabpd, ntrf, ntro3, nzjac1, & reorder, spfjsize_max, maxterms, maxfterms, & nposterms, nnegterms, nfracterms, posterms, negterms, & @@ -147,7 +147,6 @@ SUBROUTINE setup_spfuljac() INTEGER :: itemp1 INTEGER :: activity(jpcspf) -INTEGER :: permuted_nonzero_map(jpcspf,jpcspf) ! permuted map of nonzeros INTEGER, ALLOCATABLE :: permute(:,:) ! permutation matrix INTEGER, ALLOCATABLE :: map(:,:) ! a matrix of 1's indicating where the ! dense Jacobian A(i,j) has nonzero From 013a5e6d579155ccc2399eab83c9978910b5a842 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Fri, 23 Jan 2026 10:58:23 +0000 Subject: [PATCH 4/7] Add comments back in --- src/science/core/chemistry/asad/asad_sparse_vars.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/science/core/chemistry/asad/asad_sparse_vars.F90 b/src/science/core/chemistry/asad/asad_sparse_vars.F90 index 5515e39..ab2e733 100644 --- a/src/science/core/chemistry/asad/asad_sparse_vars.F90 +++ b/src/science/core/chemistry/asad/asad_sparse_vars.F90 @@ -261,6 +261,7 @@ SUBROUTINE setup_spfuljac() DO i = 1, jpcspf-1 DO j = i+1, jpcspf IF (activity(i) > activity(j)) THEN + ! exchange i and j tracers if i is more active than j. itemp1 = reorder(i) reorder(i) = reorder(j) reorder(j) = itemp1 @@ -322,6 +323,7 @@ SUBROUTINE setup_spfuljac() END DO END DO +! ! Perform error check outside of the loop to better suit GPU runs IF (total1 > spfjsize_max) THEN errcode = total1 WRITE(umMessage,'(A,2I4)') 'Total1 exceeded spfjsize_max: ', total1, & @@ -463,6 +465,7 @@ SUBROUTINE spfuljac(n_points, cdt, min_pivot, nonzero_map, spfj) INTEGER, INTENT(IN) :: n_points REAL, INTENT(IN) :: cdt REAL, INTENT(IN) :: min_pivot +INTEGER, INTENT(IN) :: nonzero_map(jpcspf,jpcspf) REAL, INTENT(OUT) :: spfj(n_points,spfjsize_max) ! Local variables From 9897a13a2c797b872556d2d7149535e52bdbc88f Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Fri, 23 Jan 2026 10:59:05 +0000 Subject: [PATCH 5/7] Add name to CONTRIBUTORS.md --- CONTRIBUTORS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index c34eb45..9e781e6 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -3,3 +3,4 @@ | GitHub user | Real Name | Affiliation | Date | | ----------- | --------- | ----------- | ---- | | james-bruten-mo | James Bruten | Met Office | 2025-12-09 | +| joewallwork | Joseph Wallwork | University of Cambridge | 2026-01-23 | From 7947249fb0cc010b71365f5ce5b2ced1c54da6d1 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Fri, 23 Jan 2026 12:38:16 +0000 Subject: [PATCH 6/7] More post-merge fixes --- src/science/core/chemistry/asad/asad_sparse_vars.F90 | 8 +++++--- src/science/core/chemistry/asad/asad_spimpmjp.F90 | 8 +++++--- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/science/core/chemistry/asad/asad_sparse_vars.F90 b/src/science/core/chemistry/asad/asad_sparse_vars.F90 index ab2e733..7185ace 100644 --- a/src/science/core/chemistry/asad/asad_sparse_vars.F90 +++ b/src/science/core/chemistry/asad/asad_sparse_vars.F90 @@ -281,6 +281,7 @@ SUBROUTINE setup_spfuljac() END DO ! Determine permutation matrix P +IF (.NOT. ALLOCATED(permute)) ALLOCATE(permute(jpcspf,jpcspf)) permute(:,:) = 0 DO i = 1, jpcspf permute(i,reorder(i)) = 1 @@ -323,7 +324,7 @@ SUBROUTINE setup_spfuljac() END DO END DO -! ! Perform error check outside of the loop to better suit GPU runs +! Perform error check outside of the loop to better suit GPU runs IF (total1 > spfjsize_max) THEN errcode = total1 WRITE(umMessage,'(A,2I4)') 'Total1 exceeded spfjsize_max: ', total1, & @@ -455,7 +456,7 @@ SUBROUTINE spfuljac(n_points, cdt, min_pivot, nonzero_map, spfj) jpspec, linfam, madvtr, moffam, ndepd, ndepw, njcoth, & nmsjac, nodd, nsjac1, nstst, ntro3, prk, spfjsize_max, & nposterms, nnegterms, nfracterms, posterms, negterms, & - fracterms, base_tracer, ffrac, y, total, nonzero_map + fracterms, base_tracer, ffrac, y, total USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook @@ -757,7 +758,7 @@ SUBROUTINE spresolv2(n_points, bb, xx, min_pivot, modified_map, spfj, max_val) ! Part (c) back-substitution: find z where U z = w ! Part (d) determine x, apply transpose(P) to z, P'z = x -USE asad_mod, ONLY: jpcspf, modified_map, reorder, spfjsize_max +USE asad_mod, ONLY: jpcspf, reorder, spfjsize_max USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook @@ -768,6 +769,7 @@ SUBROUTINE spresolv2(n_points, bb, xx, min_pivot, modified_map, spfj, max_val) REAL, INTENT(IN) :: bb(n_points,jpcspf) REAL, INTENT(OUT) :: xx(n_points,jpcspf) REAL, INTENT(IN) :: min_pivot +INTEGER, INTENT(IN) :: modified_map(jpcspf,jpcspf) REAL, INTENT(IN) :: spfj(n_points,spfjsize_max) ! Maximum tolerated value for use in filtering step. Unused if negative diff --git a/src/science/core/chemistry/asad/asad_spimpmjp.F90 b/src/science/core/chemistry/asad/asad_spimpmjp.F90 index 46d7523..86e46ff 100644 --- a/src/science/core/chemistry/asad/asad_spimpmjp.F90 +++ b/src/science/core/chemistry/asad/asad_spimpmjp.F90 @@ -256,7 +256,8 @@ SUBROUTINE asad_spimpmjp(exit_code, ix, jy, nlev, n_points, location, & solver_iter) USE asad_mod, ONLY: ptol, peps, cdt, f, fdot, nitnr, nstst, y, & - fj, nonzero_map, ltrig, jpcspf, spfj + fj, nonzero_map, ltrig, jpcspf, spfj, & + modified_map, permuted_nonzero_map USE asad_sparse_vars, ONLY: setup_spfuljac, spfuljac, spresolv2, splinslv2 USE ukca_config_specification_mod, ONLY: ukca_config USE yomhook, ONLY: lhook, dr_hook @@ -375,7 +376,7 @@ SUBROUTINE asad_spimpmjp(exit_code, ix, jy, nlev, n_points, location, & END IF !$OMP END CRITICAL (setup_jacobian_init) -CALL spfuljac(n_points,cdt,f_min,spfj) +CALL spfuljac(n_points,cdt,f_min,nonzero_map,spfj) ! Call forward Euler to make first guess, f_0 for f(t=n+1) (next timestep) CALL forward_euler(n_points, f, f_initial, f_min, nonzero_map, spfj) @@ -469,7 +470,8 @@ SUBROUTINE asad_spimpmjp(exit_code, ix, jy, nlev, n_points, location, & END DO END IF - CALL splinslv2(n_points,G_f,f_incr,f_min,f_max,spfj) + CALL splinslv2(n_points,G_f,f_incr,f_min,f_max,permuted_nonzero_map, & + modified_map,spfj) IF (ltrig .AND. printstatus == PrStatus_Diag) THEN DO jl=1,n_points From 0438f0cc4684c4f7d4be48cfbc4e733b37bf5fb9 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Fri, 23 Jan 2026 12:40:17 +0000 Subject: [PATCH 7/7] More post-merge fixes --- src/science/core/chemistry/asad/asad_sparse_vars.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/science/core/chemistry/asad/asad_sparse_vars.F90 b/src/science/core/chemistry/asad/asad_sparse_vars.F90 index 7185ace..a48a9f7 100644 --- a/src/science/core/chemistry/asad/asad_sparse_vars.F90 +++ b/src/science/core/chemistry/asad/asad_sparse_vars.F90 @@ -628,6 +628,8 @@ SUBROUTINE splinslv2(n_points, bb, xx, min_pivot, max_val, & REAL, INTENT(OUT) :: xx(n_points,jpcspf) REAL, INTENT(IN) :: min_pivot REAL, INTENT(IN) :: max_val +INTEGER, INTENT(IN) :: permuted_nonzero_map(jpcspf,jpcspf) +INTEGER, INTENT(IN OUT) :: modified_map(jpcspf,jpcspf) REAL, INTENT(IN OUT) :: spfj(1:n_points,1:spfjsize_max) ! Local variables