Changeset 4143 for LMDZ6/trunk/libf/dyn3dmem/check_isotopes_loc.F90

Timestamp:

May 9, 2022, 12:35:40 PM (2 years ago)

Author:

dcugnet

Message:

• Some variables are renamed or replaced by direct equivalents:
  • iso_indnum -> tracers(:)%iso_iName
  • niso_possibles -> niso
  • iqiso -> iqIsoPha ; index_trac -> itZonIso
  • ok_iso_verif -> isoCheck
  • ntraceurs_zone -> nzone ; ntraciso -> ntiso
  • qperemin -> min_qparent ; masseqmin -> min_qmass ; ratiomin -> min_ratio
• Some renamed variables are only aliased with the older name (using USE <module>, ONLY: <oldName> => <newName>) in routines where they are repeated many times.
• Few hard-coded indexes are now computed (examples: ilic, iso, ivap, irneb, iq_vap, iq_liq, iso_H2O, iso_HDO, iso_HTO, iso_O17, iso_O18).
• The IF(isoCheck) test is now embedded in the check_isotopes_seq and check_isotopes_loc routines (lighter calling).

File:

• LMDZ6/trunk/libf/dyn3dmem/check_isotopes_loc.F90 (moved) (moved from LMDZ6/trunk/libf/dyn3dmem/check_isotopes_loc.F) (1 diff)

LMDZ6/trunk/libf/dyn3dmem/check_isotopes_loc.F90

@@ -1 @@
-        subroutine check_isotopes(q,ijb,ije,err_msg)
-        USE infotrac, ONLY: nqtot, nqo, niso, ntraciso, nzone,
-     &                use_iso, ntraceurs_zone,
-     &                iqiso, indnum_fn_num, index_trac, tnat
-        USE parallel_lmdz
-        implicit none
+SUBROUTINE check_isotopes(q, ijb, ije, err_msg)
+   USE parallel_lmdz
+   USE strings_mod, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str
+   USE infotrac,    ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, &
+                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso, tnat
+   IMPLICIT NONE
+   include "dimensions.h"
+   REAL,             INTENT(INOUT) :: q(ijb_u:ije_u,llm,nqtot)
+   INTEGER,          INTENT(IN)    :: ijb, ije   !--- Can be local and different from ijb_u,ije_u, for example in qminimum
+   CHARACTER(LEN=*), INTENT(IN)    :: err_msg    !--- Error message to display
+   CHARACTER(LEN=maxlen) :: modname, msg1, nm(2)
+   INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar
+   INTEGER, ALLOCATABLE :: ix(:)
+   REAL    :: xtractot, xiiso, deltaD, q1, q2
+   REAL, PARAMETER :: borne     = 1e19,  &
+                      errmax    = 1e-8,  &       !--- Max. absolute error
+                      errmaxrel = 1e-3,  &       !--- Max. relative error
+                      qmin      = 1e-11, &
+                      deltaDmax =1000.0, &
+                      deltaDmin =-999.0, &
+                      ridicule  = 1e-12
+   INTEGER, SAVE :: iso_eau, iso_HDO, iso_O18, & !--- OpenMP shared variables
+                             iso_O17, iso_HTO
+   LOGICAL, SAVE :: first=.TRUE.
+!$OMP THREADPRIVATE(first)
 
-#include "dimensions.h"
+   modname='check_isotopes'
+   IF(.NOT.isoCheck)    RETURN                   !--- No need to check => finished
+   IF(isoSelect('H2O')) RETURN                   !--- No H2O isotopes group found
+   IF(niso == 0)        RETURN                   !--- No isotopes => finished
+   IF(first) THEN
+!$OMP MASTER
+      iso_eau = strIdx(isoName,'H2[16]O')
+      iso_HDO = strIdx(isoName,'H[2]HO')
+      iso_O18 = strIdx(isoName,'H2[18]O')
+      iso_O17 = strIdx(isoName,'H2[17]O')
+      iso_HTO = strIdx(isoName,'H[3]HO')
+!$OMP END MASTER
+!$OMP BARRIER
+      first = .FALSE.
+   END IF
+   CALL msg('31: err_msg='//TRIM(err_msg), modname)
 
-        ! inputs
-        integer ijb,ije ! peut être local et différent de ijb_u,ije_u, ex: dans qminimum
-        real q(ijb_u:ije_u,llm,nqtot)
-        character*(*) err_msg ! message d''erreur à afficher
+   !--- CHECK FOR NaNs (FOR ALL THE ISOTOPES, INCLUDING GEOGRAPHIC TAGGING TRACERS)
+   modname = 'check_isotopes:iso_verif_noNaN'
+   DO ixt = 1, ntiso
+      DO ipha = 1, nphas
+         iq = iqIsoPha(ixt,ipha)
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+         DO k = 1, llm
+            DO i = ijb, ije
+               IF(ABS(q(i,k,iq))<=borne) CYCLE
+               WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq)
+               CALL msg(msg1, modname)
+               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
+            END DO
+         END DO
+!$OMP END DO NOWAIT
+      END DO
+   END DO
 
-        ! locals
-        integer ixt,phase,k,i,iq,iiso,izone,ieau,iqeau
-        real xtractot,xiiso
-        real borne
-        real qmin
-        real errmax ! erreur maximale en absolu.
-        real errmaxrel ! erreur maximale en relatif autorisée
-        real deltaDmax,deltaDmin
-        real ridicule
-        parameter (borne=1e19)
-        parameter (errmax=1e-8)
-        parameter (errmaxrel=1e-3)
-        parameter (qmin=1e-11)
-        parameter (deltaDmax=1000.0,deltaDmin=-999.0)
-        parameter (ridicule=1e-12)
-        real deltaD
+   !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL)
+   modname = 'check_isotopes:iso_verif_egalite'
+   ixt = iso_eau
+   IF(ixt /= 0) THEN
+      DO ipha = 1, nphas
+         iq = iqIsoPha(ixt,ipha)
+         iqpar = tracers(iq)%iqParent
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+         DO k = 1, llm
+            DO i = ijb, ije
+               q1 = q(i,k,iqpar)
+               q2 = q(i,k,iq)
+!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
+!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
+!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
+!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
+               IF(ABS(q1-q2) <= errmax .OR. ABS(q1-q2)/MAX(MAX(ABS(q1),ABS(q2)),1e-18) <= errmaxrel) THEN
+                  q(i,k,iq) = q1                 !--- Bidouille pour convergence
+!                 q(i,k,tracers(iqPar)%iqDesc) = q(i,k,tracers(iqPar)%iqDesc) * q1 / q2
+                  CYCLE
+               END IF
+               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
+               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
+               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
+               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
+               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
+            END DO
+         END DO
+!$OMP END DO NOWAIT
+      END DO
+   END IF
 
-        if (niso > 0) then
+   !--- CHECK DELTA ANOMALIES
+   modname = 'check_isotopes:iso_verif_aberrant'
+   ix = [ iso_HDO  ,   iso_O18 ]
+   nm = ['deltaD  ', 'deltaO18']
+   DO iiso = 1, SIZE(ix)
+      ixt = ix(iiso)
+      IF(ixt  == 0) CYCLE
+      DO ipha = 1, nphas
+         iq = iqIsoPha(ixt,ipha)
+         iqpar = tracers(iq)%iqParent
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+         DO k = 1, llm
+            DO i = ijb, ije
+               q1 = q(i,k,iqpar)
+               q2 = q(i,k,iq)
+!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
+!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
+!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
+!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
+               IF(q2 <= qmin) CYCLE
+               deltaD = (q2/q1/tnat(ixt)-1.)*1000.
+               IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
+               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
+               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
+               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
+               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
+               CALL msg(TRIM(nm(iiso))//TRIM(real2str(deltaD)), modname)
+               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
+            END DO
+         END DO
+!$OMP END DO NOWAIT
+      END DO
+   END DO
 
-!        write(*,*) 'check_isotopes 31: err_msg=',err_msg
-        ! verifier que rien n'est NaN
-        do ixt=1,ntraciso
-          do phase=1,nqo
-            iq=iqiso(ixt,phase)
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-            do k=1,llm
-              DO i = ijb,ije
-                if ((q(i,k,iq).gt.-borne).and.
-     :            (q(i,k,iq).lt.borne)) then
-                else !if ((x(ixt,i,j).gt.-borne).and.
-                  write(*,*) 'erreur detectee par iso_verif_noNaN:'
-                  write(*,*) err_msg
-                  write(*,*) 'q,i,k,iq=',q(i,k,iq),i,k,iq
-                  write(*,*) 'borne=',borne
-                  call abort_gcm('check_isotopes_loc','plantage iso',0)
-                endif  !if ((x(ixt,i,j).gt.-borne).and.
-              enddo !DO i = ijb,ije
-            enddo !do k=1,llm
-c$OMP END DO NOWAIT
-          enddo !do phase=1,nqo
-        enddo !do ixt=1,ntraciso
+   IF(nzone == 0) RETURN
 
-!        write(*,*) 'check_isotopes 52'
-        ! verifier que l'eau normale est OK
-        if (use_iso(1)) then
-          ixt=indnum_fn_num(1)
-          do phase=1,nqo
-            iq=iqiso(ixt,phase)
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-            do k=1,llm
-            DO i = ijb,ije  
-              if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
-     :          (abs((q(i,k,phase)-q(i,k,iq))/
-     :           max(max(abs(q(i,k,phase)),abs(q(i,k,iq))),1e-18))
-     :           .gt.errmaxrel)) then
-                  write(*,*) 'erreur detectee par iso_verif_egalite:'
-                  write(*,*) err_msg
-                  write(*,*) 'ixt,phase,ijb=',ixt,phase,ijb
-                  write(*,*) 'q,iq,i,k=',q(i,k,iq),iq,i,k
-                  write(*,*) 'q(i,k,phase)=',q(i,k,phase)
-                  call abort_gcm('check_isotopes_loc','plantage iso',0)
-              endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
-              ! bidouille pour éviter divergence:
-              q(i,k,iq)= q(i,k,phase) 
-            enddo ! DO i = ijb,ije
-            enddo !do k=1,llm
-c$OMP END DO NOWAIT
-          enddo ! do phase=1,nqo 
-        endif !if (use_iso(1)) then
-        
-!        write(*,*) 'check_isotopes 78'
-        ! verifier que HDO est raisonable
-        if (use_iso(2)) then
-          ixt=indnum_fn_num(2)
-          do phase=1,nqo
-            iq=iqiso(ixt,phase)
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-            do k=1,llm
-            DO i = ijb,ije
-            if (q(i,k,iq).gt.qmin) then
-             deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(2)-1)*1000
-             if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
-                  write(*,*) 'erreur detectee par iso_verif_aberrant:'
-                  write(*,*) err_msg
-                  write(*,*) 'ixt,phase=',ixt,phase
-                  write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k
-                  write(*,*) 'q=',q(i,k,:)
-                  write(*,*) 'deltaD=',deltaD
-                  call abort_gcm('check_isotopes_loc','plantage iso',0)
-             endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
-            endif !if (q(i,k,iq).gt.qmin) then
-            enddo !DO i = ijb,ije
-            enddo !do k=1,llm
-c$OMP END DO NOWAIT
-          enddo ! do phase=1,nqo 
-        endif !if (use_iso(2)) then
+   !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES
+   modname = 'check_isotopes:iso_verif_aberrant'
+   IF(iso_eau /= 0 .AND. iso_HDO /= 0) THEN
+      DO izon = 1, nzone
+         ixt  = itZonIso(izon, iso_HDO)
+         ieau = itZonIso(izon, iso_eau)
+         DO ipha = 1, nphas
+            iq    = iqIsoPha(ixt,  ipha)
+            iqeau = iqIsoPha(ieau, ipha)
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+            DO k = 1, llm
+               DO i = ijb, ije
+                  q1 = q(i,k,iqeau)
+                  q2 = q(i,k,iq)
+                  IF(q2<=qmin) CYCLE
+                  deltaD = (q2/q1/tnat(iso_HDO)-1.)*1000.
+                  IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
+                  CALL msg('izon, ipha = '//TRIM(strStack(int2str([izon, ipha]))), modname)
+                  CALL msg( 'ixt, ieau = '//TRIM(strStack(int2str([ ixt, ieau]))), modname)
+                  msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
+                  CALL msg(TRIM(tracers(iqeau)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
+                  CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
+                  CALL msg('deltaD = '//TRIM(real2str(deltaD)), modname)
+                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
+               END DO
+            END DO
+!$OMP END DO NOWAIT
+         END DO
+      END DO
+   END IF
 
-!        write(*,*) 'check_isotopes 103'
-        ! verifier que O18 est raisonable
-        if (use_iso(3)) then
-          ixt=indnum_fn_num(3)
-          do phase=1,nqo
-            iq=iqiso(ixt,phase)
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-            do k=1,llm
-            DO i = ijb,ije
-            if (q(i,k,iq).gt.qmin) then
-             deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(3)-1)*1000
-             if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
-                  write(*,*) 'erreur detectee iso_verif_aberrant O18:'
-                  write(*,*) err_msg
-                  write(*,*) 'ixt,phase=',ixt,phase
-                  write(*,*) 'q,iq,i,k,=',q(i,k,phase),iq,i,k
-                  write(*,*) 'xt=',q(i,k,:)
-                  write(*,*) 'deltaO18=',deltaD
-                  call abort_gcm('check_isotopes_loc','plantage iso',0)
-             endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
-            endif !if (q(i,k,iq).gt.qmin) then
-            enddo !DO i = ijb,ije
-            enddo !do k=1,llm
-c$OMP END DO NOWAIT
-          enddo ! do phase=1,nqo 
-        endif !if (use_iso(2)) then
+   !--- CHECK FOR TAGGING TRACERS CONSERVATION (PARENT AND TAGGING TRACERS SUM OVER ALL REGIONS MUST BE EQUAL)
+   DO iiso = 1, niso
+      DO ipha = 1, nphas
+         iq = iqIsoPha(iiso, ipha)
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+         DO k = 1, llm
+            DO i = ijb, ije
+               xiiso = q(i,k,iq)
+               xtractot = SUM(q(i, k, iqIsoPha(itZonIso(1:nzone,iiso), ipha)))
+               IF(ABS(xtractot-xiiso) > errmax .AND. ABS(xtractot-xiiso)/MAX(MAX(ABS(xtractot),ABS(xiiso)),1e-18) > errmaxrel) THEN
+                  CALL msg('Error in iso_verif_aberrant trac: '//TRIM(err_msg))
+                  CALL msg('iiso, ipha = '//TRIM(strStack(int2str([iiso, ipha]))), modname)
+                  CALL msg('q('//TRIM(strStack(int2str([i,k])))//',:) = '//TRIM(strStack(real2str(q(i,k,:)))), modname)
+                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
+               END IF
+               IF(ABS(xtractot) <= ridicule) CYCLE
+               DO izon = 1, nzone
+                  q(i,k,iq) = q(i,k,iq) / xtractot * xiiso !--- Bidouille pour convergence
+               END DO
+            END DO
+         END DO
+!$OMP END DO NOWAIT
+      END DO
+   END DO
 
+END SUBROUTINE check_isotopes
 
-!        write(*,*) 'check_isotopes 129'
-        if (nzone > 0) then
-
-          if (use_iso(2).and.use_iso(1)) then
-            do izone=1,ntraceurs_zone
-             ixt=index_trac(izone,indnum_fn_num(2))
-             ieau=index_trac(izone,indnum_fn_num(1))
-             do phase=1,nqo
-               iq=iqiso(ixt,phase)
-               iqeau=iqiso(ieau,phase)
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-               do k=1,llm
-                DO i = ijb,ije
-                if (q(i,k,iq).gt.qmin) then
-                 deltaD=(q(i,k,iq)/q(i,k,iqeau)/tnat(2)-1)*1000
-                 if ((deltaD.gt.deltaDmax).or.
-     &                   (deltaD.lt.deltaDmin)) then
-                  write(*,*) 'erreur dans iso_verif_aberrant trac:'
-                  write(*,*) err_msg
-                  write(*,*) 'izone,phase=',izone,phase
-                  write(*,*) 'ixt,ieau=',ixt,ieau
-                  write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k
-                  write(*,*) 'deltaD=',deltaD                  
-                  call abort_gcm('check_isotopes_loc','plantage iso',0)
-                 endif !if ((deltaD.gt.deltaDmax).or.
-                endif !if (q(i,k,iq).gt.qmin) then
-                enddo !DO i = ijb,ije
-                enddo  ! do k=1,llm
-c$OMP END DO NOWAIT
-              enddo ! do phase=1,nqo    
-            enddo !do izone=1,ntraceurs_zone
-          endif !if (use_iso(2).and.use_iso(1)) then
-
-          do iiso=1,niso
-           do phase=1,nqo
-              iq=iqiso(iiso,phase)
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-              do k=1,llm
-                DO i = ijb,ije
-                   xtractot=0.0
-                   xiiso=q(i,k,iq)
-                   do izone=1,ntraceurs_zone
-                      iq=iqiso(index_trac(izone,iiso),phase)
-                      xtractot=xtractot+ q(i,k,iq)
-                   enddo !do izone=1,ntraceurs_zone
-                   if ((abs(xtractot-xiiso).gt.errmax).and.
-     :                  (abs(xtractot-xiiso)/
-     :                  max(max(abs(xtractot),abs(xiiso)),1e-18)
-     :                  .gt.errmaxrel)) then
-                  write(*,*) 'erreur detectee par iso_verif_traceurs:'
-                  write(*,*) err_msg
-                  write(*,*) 'iiso,phase=',iiso,phase
-                  write(*,*) 'i,k,=',i,k
-                  write(*,*) 'q(i,k,:)=',q(i,k,:)
-                  call abort_gcm('check_isotopes_loc','plantage iso',0)
-                 endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
-                  
-                 ! bidouille pour éviter divergence:
-                 if (abs(xtractot).gt.ridicule) then
-                   do izone=1,ntraceurs_zone
-                     ixt=index_trac(izone,iiso) 
-                     q(i,k,iq)=q(i,k,iq)/xtractot*xiiso
-                   enddo !do izone=1,ntraceurs_zone                
-                  endif !if ((abs(xtractot).gt.ridicule) then
-                enddo !DO i = ijb,ije
-              enddo !do k=1,llm
-c$OMP END DO NOWAIT
-           enddo !do phase=1,nqo
-          enddo !do iiso=1,niso
-
-        endif !if (nzone > 0)
-
-        endif ! if (niso > 0)
-!        write(*,*) 'check_isotopes 198'
-        
-        end
-
-