Ignore:
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:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/dyn3d/check_isotopes.F90

    r4142 r4143  
    1         subroutine check_isotopes_seq(q,ip1jmp1,err_msg)
    2         USE infotrac, ONLY: nqtot, nqo, niso, ntraciso, nzone,
    3      &                     use_iso,
    4      &                     iqiso, index_trac,indnum_fn_num, tnat
    5         implicit none
     1SUBROUTINE check_isotopes_seq(q, ip1jmp1, err_msg)
     2   USE strings_mod, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str
     3   USE infotrac,    ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, &
     4                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso, tnat
     5   IMPLICIT NONE
     6   include "dimensions.h"
     7   REAL,             INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
     8   INTEGER,          INTENT(IN)    :: ip1jmp1
     9   CHARACTER(LEN=*), INTENT(IN)    :: err_msg    !--- Error message to display
     10   CHARACTER(LEN=maxlen) :: modname, msg1, nm(2)
     11   INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar
     12   INTEGER, ALLOCATABLE :: ix(:)
     13   REAL    :: xtractot, xiiso, deltaD, q1, q2
     14   REAL, PARAMETER :: borne     = 1e19,  &
     15                      errmax    = 1e-8,  &       !--- Max. absolute error
     16                      errmaxrel = 1e-3,  &       !--- Max. relative error
     17                      qmin      = 1e-11, &
     18                      deltaDmax =1000.0, &
     19                      deltaDmin =-999.0, &
     20                      ridicule  = 1e-12
     21   INTEGER, SAVE :: iso_eau, iso_HDO, iso_O18, &
     22                             iso_O17, iso_HTO
     23   LOGICAL, SAVE :: first=.TRUE.
    624
    7 #include "dimensions.h"
     25   modname='check_isotopes'
     26   IF(.NOT.isoCheck)    RETURN                   !--- No need to check => finished
     27   IF(isoSelect('H2O')) RETURN                   !--- No H2O isotopes group found
     28   IF(niso == 0)        RETURN                   !--- No isotopes => finished
     29   IF(first) THEN
     30      iso_eau = strIdx(isoName,'H2[16]O')
     31      iso_HDO = strIdx(isoName,'H[2]HO')
     32      iso_O18 = strIdx(isoName,'H2[18]O')
     33      iso_O17 = strIdx(isoName,'H2[17]O')
     34      iso_HTO = strIdx(isoName,'H[3]HO')
     35      first = .FALSE.
     36   END IF
     37   CALL msg('31: err_msg='//TRIM(err_msg), modname)
    838
    9         ! inputs
    10         integer ip1jmp1
    11         real q(ip1jmp1,llm,nqtot)
    12         character*(*) err_msg ! message d''erreur à afficher
     39   !--- CHECK FOR NaNs (FOR ALL THE ISOTOPES, INCLUDING GEOGRAPHIC TAGGING TRACERS)
     40   modname = 'check_isotopes:iso_verif_noNaN'
     41   DO ixt = 1, ntiso
     42      DO ipha = 1, nphas
     43         iq = iqIsoPha(ixt,ipha)
     44         DO k = 1, llm
     45            DO i = 1, ip1jmp1
     46               IF(ABS(q(i,k,iq)) < borne) CYCLE
     47               WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq)
     48               CALL msg(msg1, modname)
     49               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
     50            END DO
     51         END DO
     52      END DO
     53   END DO
    1354
    14         ! locals
    15         integer ixt,phase,k,i,iq,iiso,izone,ieau,iqeau
    16         real xtractot,xiiso
    17         real borne
    18         real qmin
    19         real errmax ! erreur maximale en absolu.
    20         real errmaxrel ! erreur maximale en relatif autorisée
    21         real deltaDmax,deltaDmin
    22         real ridicule
    23         parameter (borne=1e19)
    24         parameter (errmax=1e-8)
    25         parameter (errmaxrel=1e-3)
    26         parameter (qmin=1e-11)
    27         parameter (deltaDmax=1000.0,deltaDmin=-999.0)
    28         parameter (ridicule=1e-12)
    29         real deltaD
     55   !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL)
     56   modname = 'check_isotopes:iso_verif_egalite'
     57   ixt = iso_eau
     58   IF(ixt /= 0) THEN
     59      DO ipha = 1, nphas
     60         iq = iqIsoPha(ixt,ipha)
     61         iqpar = tracers(iq)%iqParent
     62         DO k = 1, llm
     63            DO i = 1, ip1jmp1
     64               q1 = q(i,k,iqpar)
     65               q2 = q(i,k,iq)
     66!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
     67!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
     68!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
     69!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
     70               IF(ABS(q1-q2) <= errmax .OR. ABS(q1-q2)/MAX(MAX(ABS(q1),ABS(q2)),1e-18) <= errmaxrel) THEN
     71                  q(i,k,iq) = q1                 !--- Bidouille pour convergence
     72!                 q(i,k,tracers(iqPar)%iqDesc) = q(i,k,tracers(iqPar)%iqDesc) * q1 / q2
     73                  CYCLE
     74               END IF
     75               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
     76               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
     77               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
     78               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
     79               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
     80            END DO
     81         END DO
     82      END DO
     83   END IF
    3084
    31         if (niso > 0) then
     85   !--- CHECK DELTA ANOMALIES
     86   modname = 'check_isotopes:iso_verif_aberrant'
     87   ix = [ iso_HDO  ,   iso_O18 ]
     88   nm = ['deltaD  ', 'deltaO18']
     89   DO iiso = 1, SIZE(ix)
     90      ixt = ix(iiso)
     91      IF(ixt  == 0) CYCLE
     92      DO ipha = 1, nphas
     93         iq = iqIsoPha(ixt,ipha)
     94         iqpar = tracers(iq)%iqParent
     95         DO k = 1, llm
     96            DO i = 1, ip1jmp1
     97               q1 = q(i,k,iqpar)
     98               q2 = q(i,k,iq)
     99!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
     100!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
     101!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
     102!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
     103               IF(q2 <= qmin) CYCLE
     104               deltaD = (q2/q1/tnat(ixt)-1.)*1000.
     105               IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
     106               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
     107               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
     108               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
     109               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
     110               CALL msg(TRIM(nm(iiso))//TRIM(real2str(deltaD)), modname)
     111               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
     112            END DO
     113         END DO
     114      END DO
     115   END DO
    32116
    33         write(*,*) 'check_isotopes 31: err_msg=',err_msg
    34         ! verifier que rien n'est NaN
    35         do ixt=1,ntraciso
    36           do phase=1,nqo
    37             iq=iqiso(ixt,phase)
    38             do k=1,llm
    39               DO i = 1,ip1jmp1
    40                 if ((q(i,k,iq).gt.-borne).and.
    41      :            (q(i,k,iq).lt.borne)) then
    42                 else !if ((x(ixt,i,j).gt.-borne).and.
    43                   write(*,*) 'erreur detectee par iso_verif_noNaN:'
    44                   write(*,*) err_msg
    45                   write(*,*) 'q,i,k,iq=',q(i,k,iq),i,k,iq
    46                   write(*,*) 'borne=',borne
    47                   call abort_gcm('check_isotopes_loc','plantage iso',0)
    48                 endif  !if ((x(ixt,i,j).gt.-borne).and.
    49               enddo !DO i = 1,ip1jmp1
    50             enddo !do k=1,llm
    51           enddo !do phase=1,nqo
    52         enddo !do ixt=1,ntraciso
     117   IF(nzone == 0) RETURN
    53118
    54         !write(*,*) 'check_isotopes 52'
    55         ! verifier que l'eau normale est OK
    56         if (use_iso(1)) then
    57           ixt=indnum_fn_num(1)
    58           do phase=1,nqo
    59             iq=iqiso(ixt,phase)
    60             do k=1,llm
    61             DO i = 1,ip1jmp1 
    62               if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
    63      :          (abs((q(i,k,phase)-q(i,k,iq))/
    64      :           max(max(abs(q(i,k,phase)),abs(q(i,k,iq))),1e-18))
    65      :           .gt.errmaxrel)) then
    66                   write(*,*) 'erreur detectee par iso_verif_egalite:'
    67                   write(*,*) err_msg
    68                   write(*,*) 'ixt,phase=',ixt,phase
    69                   write(*,*) 'q,iq,i,k=',q(i,k,iq),iq,i,k
    70                   write(*,*) 'q(i,k,phase)=',q(i,k,phase)
    71                   call abort_gcm('check_isotopes_loc','plantage iso',0)
    72               endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
    73               ! bidouille pour éviter divergence:
    74               q(i,k,iq)= q(i,k,phase)
    75             enddo ! DO i = 1,ip1jmp1
    76             enddo !do k=1,llm
    77           enddo ! do phase=1,nqo
    78         endif !if (use_iso(1)) then
    79        
    80         !write(*,*) 'check_isotopes 78'
    81         ! verifier que HDO est raisonable
    82         if (use_iso(2)) then
    83           ixt=indnum_fn_num(2)
    84           do phase=1,nqo
    85             iq=iqiso(ixt,phase)
    86             do k=1,llm
    87             DO i = 1,ip1jmp1
    88             if (q(i,k,iq).gt.qmin) then
    89              deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(2)-1)*1000
    90              if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
    91                   write(*,*) 'erreur detectee par iso_verif_aberrant:'
    92                   write(*,*) err_msg
    93                   write(*,*) 'ixt,phase=',ixt,phase
    94                   write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k
    95                   write(*,*) 'q=',q(i,k,:)
    96                   write(*,*) 'deltaD=',deltaD
    97                   call abort_gcm('check_isotopes_loc','plantage iso',0)
    98              endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
    99             endif !if (q(i,k,iq).gt.qmin) then
    100             enddo !DO i = 1,ip1jmp1
    101             enddo !do k=1,llm
    102           enddo ! do phase=1,nqo
    103         endif !if (use_iso(2)) then
     119   !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES
     120   modname = 'check_isotopes:iso_verif_aberrant'
     121   IF(iso_eau /= 0 .AND. iso_HDO /= 0) THEN
     122      DO izon = 1, nzone
     123         ixt  = itZonIso(izon, iso_HDO)
     124         ieau = itZonIso(izon, iso_eau)
     125         DO ipha = 1, nphas
     126            iq    = iqIsoPha(ixt,  ipha)
     127            iqeau = iqIsoPha(ieau, ipha)
     128            DO k = 1, llm
     129               DO i = 1, ip1jmp1
     130                  q1 = q(i,k,iqeau)
     131                  q2 = q(i,k,iq)
     132                  IF(q2<=qmin) CYCLE
     133                  deltaD = (q2/q1/tnat(iso_HDO)-1.)*1000.
     134                  IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
     135                  CALL msg('izon, ipha = '//TRIM(strStack(int2str([izon, ipha]))), modname)
     136                  CALL msg( 'ixt, ieau = '//TRIM(strStack(int2str([ ixt, ieau]))), modname)
     137                  msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
     138                  CALL msg(TRIM(tracers(iqeau)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
     139                  CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
     140                  CALL msg('deltaD = '//TRIM(real2str(deltaD)), modname)
     141                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
     142               END DO
     143            END DO
     144         END DO
     145      END DO
     146   END IF
    104147
    105         !write(*,*) 'check_isotopes 103'
    106         ! verifier que O18 est raisonable
    107         if (use_iso(3)) then
    108           ixt=indnum_fn_num(3)
    109           do phase=1,nqo
    110             iq=iqiso(ixt,phase)
    111             do k=1,llm
    112             DO i = 1,ip1jmp1
    113             if (q(i,k,iq).gt.qmin) then
    114              deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(3)-1)*1000
    115              if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
    116                   write(*,*) 'erreur detectee iso_verif_aberrant O18:'
    117                   write(*,*) err_msg
    118                   write(*,*) 'ixt,phase=',ixt,phase
    119                   write(*,*) 'q,iq,i,k,=',q(i,k,phase),iq,i,k
    120                   write(*,*) 'xt=',q(i,k,:)
    121                   write(*,*) 'deltaO18=',deltaD
    122                   call abort_gcm('check_isotopes_loc','plantage iso',0)
    123              endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
    124             endif !if (q(i,k,iq).gt.qmin) then
    125             enddo !DO i = 1,ip1jmp1
    126             enddo !do k=1,llm
    127           enddo ! do phase=1,nqo
    128         endif !if (use_iso(2)) then
     148   !--- CHECK FOR TAGGING TRACERS CONSERVATION (PARENT AND TAGGING TRACERS SUM OVER ALL REGIONS MUST BE EQUAL)
     149   DO iiso = 1, niso
     150      DO ipha = 1, nphas
     151         iq = iqIsoPha(iiso, ipha)
     152         DO k = 1, llm
     153            DO i = 1, ip1jmp1
     154               xiiso = q(i,k,iq)
     155               xtractot = SUM(q(i, k, iqIsoPha(itZonIso(1:nzone,iiso), ipha)))
     156               IF(ABS(xtractot-xiiso) > errmax .AND. ABS(xtractot-xiiso)/MAX(MAX(ABS(xtractot),ABS(xiiso)),1e-18) > errmaxrel) THEN
     157                  CALL msg('Error in iso_verif_aberrant trac: '//TRIM(err_msg))
     158                  CALL msg('iiso, ipha = '//TRIM(strStack(int2str([iiso, ipha]))), modname)
     159                  CALL msg('q('//TRIM(strStack(int2str([i,k])))//',:) = '//TRIM(strStack(real2str(q(i,k,:)))), modname)
     160                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
     161               END IF
     162               IF(ABS(xtractot) <= ridicule) CYCLE
     163               DO izon = 1, nzone
     164                  q(i,k,iq) = q(i,k,iq) / xtractot * xiiso !--- Bidouille pour convergence
     165               END DO
     166            END DO
     167         END DO
     168      END DO
     169   END DO
    129170
     171END SUBROUTINE check_isotopes_seq
    130172
    131         !write(*,*) 'check_isotopes 129'
    132         if (nzone > 0) then
    133 
    134           if (use_iso(2).and.use_iso(1)) then
    135             do izone=1,nzone
    136              ixt=index_trac(izone,indnum_fn_num(2))
    137              ieau=index_trac(izone,indnum_fn_num(1))
    138              do phase=1,nqo
    139                iq=iqiso(ixt,phase)
    140                iqeau=iqiso(ieau,phase)
    141                do k=1,llm
    142                 DO i = 1,ip1jmp1
    143                 if (q(i,k,iq).gt.qmin) then
    144                  deltaD=(q(i,k,iq)/q(i,k,iqeau)/tnat(2)-1)*1000
    145                  if ((deltaD.gt.deltaDmax).or.
    146      &                   (deltaD.lt.deltaDmin)) then
    147                   write(*,*) 'erreur dans iso_verif_aberrant trac:'
    148                   write(*,*) err_msg
    149                   write(*,*) 'izone,phase=',izone,phase
    150                   write(*,*) 'ixt,ieau=',ixt,ieau
    151                   write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k
    152                   write(*,*) 'deltaD=',deltaD
    153                   call abort_gcm('check_isotopes_loc','plantage iso',0)
    154                  endif !if ((deltaD.gt.deltaDmax).or.
    155                 endif !if (q(i,k,iq).gt.qmin) then
    156                 enddo !DO i = 1,ip1jmp1
    157                 enddo  ! do k=1,llm
    158               enddo ! do phase=1,nqo   
    159             enddo !do izone=1,nzone
    160           endif !if (use_iso(2).and.use_iso(1)) then
    161 
    162           do iiso=1,niso
    163            do phase=1,nqo
    164               iq=iqiso(iiso,phase)
    165               do k=1,llm
    166                 DO i = 1,ip1jmp1
    167                    xtractot=0.0
    168                    xiiso=q(i,k,iq)
    169                    do izone=1,nzone
    170                       iq=iqiso(index_trac(izone,iiso),phase)
    171                       xtractot=xtractot+ q(i,k,iq)
    172                    enddo !do izone=1,ntraceurs_zone
    173                    if ((abs(xtractot-xiiso).gt.errmax).and.
    174      :                  (abs(xtractot-xiiso)/
    175      :                  max(max(abs(xtractot),abs(xiiso)),1e-18)
    176      :                  .gt.errmaxrel)) then
    177                   write(*,*) 'erreur detectee par iso_verif_traceurs:'
    178                   write(*,*) err_msg
    179                   write(*,*) 'iiso,phase=',iiso,phase
    180                   write(*,*) 'i,k,=',i,k
    181                   write(*,*) 'q(i,k,:)=',q(i,k,:)
    182                   call abort_gcm('check_isotopes_loc','plantage iso',0)
    183                  endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
    184                  
    185                  ! bidouille pour éviter divergence:
    186                  if (abs(xtractot).gt.ridicule) then
    187                    do izone=1,nzone
    188                      ixt=index_trac(izone,iiso)
    189                      q(i,k,iq)=q(i,k,iq)/xtractot*xiiso
    190                    enddo !do izone=1,nzone               
    191                   endif !if ((abs(xtractot).gt.ridicule) then
    192                 enddo !DO i = 1,ip1jmp1
    193               enddo !do k=1,llm
    194            enddo !do phase=1,nqo
    195           enddo !do iiso=1,niso
    196 
    197         endif !if (nzone > 0)
    198 
    199         endif ! if (niso > 0)
    200         !write(*,*) 'check_isotopes 198'
    201        
    202         end
    203 
Note: See TracChangeset for help on using the changeset viewer.