Changeset 2270 for LMDZ5/trunk/libf/dyn3d_common

Timestamp:

May 7, 2015, 5:45:04 PM (9 years ago)

Author:

crisi

Message:

Adding isotopes in the dynamics and more generally tracers of tracers.
CRisi

Location:

LMDZ5/trunk/libf/dyn3d_common

Files:

• infotrac.F90 (modified) (11 diffs)
• iso_verif_dyn.F (added)

LMDZ5/trunk/libf/dyn3d_common/infotrac.F90

@@ -12 +12 @@
   INTEGER, SAVE :: nbtr
 
+! CRisi: nb traceurs pères= directement advectés par l'air
+  INTEGER, SAVE :: nqperes
+
 ! Name variables
   CHARACTER(len=20), ALLOCATABLE, DIMENSION(:), SAVE :: tname ! tracer short name for restart and diagnostics
@@ -22 +25 @@
 !         dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code. 
   INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: niadv ! equivalent dyn / physique
+
+! CRisi: tableaux de fils
+  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: nqfils
+  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: nqdesc ! nombres de fils + nombre de tous les petits fils sur toutes les générations
+  INTEGER, SAVE :: nqdesc_tot
+  INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE    :: iqfils
+  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: iqpere
 
 ! conv_flg(it)=0 : convection desactivated for tracer number it 
@@ -30 +40 @@
   CHARACTER(len=4),SAVE :: type_trac
   CHARACTER(len=8),DIMENSION(:),ALLOCATABLE, SAVE :: solsym
+   
+    ! CRisi: cas particulier des isotopes
+    LOGICAL,SAVE :: ok_isotopes,ok_iso_verif,ok_isotrac,ok_init_iso
+    INTEGER :: niso_possibles   
+    PARAMETER ( niso_possibles=5)
+    real, DIMENSION (niso_possibles),SAVE :: tnat,alpha_ideal
+    LOGICAL, DIMENSION(niso_possibles),SAVE ::  use_iso
+    INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE ::  iqiso ! donne indice iq en fn de (ixt,phase) 
+    INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  iso_num ! donne numéro iso entre 1 et niso_possibles en fn de nqtot
+    INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  iso_indnum ! donne numéro iso entre 1 et niso effectif en fn de nqtot
+    INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  zone_num ! donne numéro de la zone de tracage en fn de nqtot
+    INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  phase_num ! donne numéro de la zone de tracage en fn de nqtot
+    INTEGER, DIMENSION(niso_possibles), SAVE :: indnum_fn_num ! donne indice entre entre 1 et niso en fonction du numéro d isotope entre 1 et niso_possibles
+    INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE ::  index_trac ! numéro ixt en fn izone, indnum entre 1 et niso
+    INTEGER,SAVE :: niso,ntraceurs_zone,ntraciso
  
 CONTAINS
@@ -63 +88 @@
 
     CHARACTER(len=15), ALLOCATABLE, DIMENSION(:) :: tnom_0  ! tracer short name
+    CHARACTER(len=15), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name: CRisi
     CHARACTER(len=3), DIMENSION(30) :: descrq
     CHARACTER(len=1), DIMENSION(3)  :: txts
@@ -70 +96 @@
     INTEGER :: nqtrue  ! number of tracers read from tracer.def, without higer order of moment
     INTEGER :: iq, new_iq, iiq, jq, ierr
+    INTEGER :: ifils,ipere,generation ! CRisi
+    LOGICAL :: continu,nouveau_traceurdef
+    INTEGER :: IOstatus ! gestion de la retrocompatibilite de traceur.def
+    CHARACTER(len=15) :: tchaine    
 
     character(len=*),parameter :: modname="infotrac_init"
@@ -134 +164 @@
           WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
           READ(90,*) nqtrue
+          write(lunout,*) 'nqtrue=',nqtrue
        ELSE 
           WRITE(lunout,*) trim(modname),': Problem in opening traceur.def'
@@ -192 +223 @@
 ! Allocate variables depending on nqtrue
 !
-    ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue))
+    ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue),tnom_transp(nqtrue))
 !
 !jyg<
@@ -230 +261 @@
           ! Continue to read tracer.def
           DO iq=1,nqtrue
-             READ(90,*) hadv(iq),vadv(iq),tnom_0(iq)
-          END DO
+
+             write(*,*) 'infotrac 237: iq=',iq
+             ! CRisi: ajout du nom du fluide transporteur
+             ! mais rester retro compatible
+             READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine
+             write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq)
+             write(lunout,*) 'tchaine=',trim(tchaine)
+             write(*,*) 'infotrac 238: IOstatus=',IOstatus
+             if (IOstatus.ne.0) then
+                CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1)
+             endif
+             ! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un
+             ! espace ou pas au milieu de la chaine.
+             continu=1
+             nouveau_traceurdef=0
+             iiq=1
+             do while (continu)
+                if (tchaine(iiq:iiq).eq.' ') then
+                  nouveau_traceurdef=1
+                  continu=0
+                else if (iiq.lt.LEN_TRIM(tchaine)) then
+                  iiq=iiq+1
+                else
+                  continu=0      
+                endif
+             enddo
+             write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef
+             if (nouveau_traceurdef) then
+                write(lunout,*) 'C''est la nouvelle version de traceur.def'
+                tnom_0(iq)=tchaine(1:iiq-1)
+                tnom_transp(iq)=tchaine(iiq+1:15)
+             else
+                write(lunout,*) 'C''est l''ancienne version de traceur.def'
+                write(lunout,*) 'On suppose que les traceurs sont tous d''air'
+                tnom_0(iq)=tchaine
+                tnom_transp(iq) = 'air'
+             endif
+             write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>'
+             write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>'
+
+          END DO !DO iq=1,nqtrue
           CLOSE(90)  
+
        ELSE ! Without tracer.def, set default values 
          if (planet_type=="earth") then
@@ -239 +310 @@
           vadv(1) = 14
           tnom_0(1) = 'H2Ov'
+          tnom_transp(1) = 'air'
           hadv(2) = 10
           vadv(2) = 10
           tnom_0(2) = 'H2Ol'
+          tnom_transp(2) = 'air'
           hadv(3) = 10
           vadv(3) = 10
           tnom_0(3) = 'RN'
+          tnom_transp(3) = 'air'
           hadv(4) = 10
           vadv(4) = 10
           tnom_0(4) = 'PB'
+          tnom_transp(4) = 'air'
          else ! default for other planets
           hadv(1) = 10
           vadv(1) = 10
           tnom_0(1) = 'dummy'
+          tnom_transp(1) = 'dummy'
          endif ! of if (planet_type=="earth")
        END IF
@@ -258 +334 @@
        WRITE(lunout,*) trim(modname),': nombre de traceurs ',nqtrue
        DO iq=1,nqtrue
-          WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq)
+          WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq),tnom_transp(iq)
        END DO
 
@@ -447 +523 @@
     END DO
 
+
+! CRisi: quels sont les traceurs fils et les traceurs pères.
+! initialiser tous les tableaux d'indices liés aux traceurs familiaux
+! + vérifier que tous les pères sont écrits en premières positions
+    ALLOCATE(nqfils(nqtot),nqdesc(nqtot))    
+    ALLOCATE(iqfils(nqtot,nqtot))   
+    ALLOCATE(iqpere(nqtot))
+    nqperes=0
+    nqfils(:)=0
+    nqdesc(:)=0
+    iqfils(:,:)=0
+    iqpere(:)=0
+    nqdesc_tot=0    
+    DO iq=1,nqtot
+      if (tnom_transp(iq) == 'air') then
+        ! ceci est un traceur père
+        WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un pere'
+        nqperes=nqperes+1
+        iqpere(iq)=0
+      else !if (tnom_transp(iq) == 'air') then
+        ! ceci est un fils. Qui est son père?
+        WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un fils'
+        continu=.true.
+        ipere=1
+        do while (continu)           
+          if (tnom_transp(iq) == tnom_0(ipere)) then
+            ! Son père est ipere
+            WRITE(lunout,*) 'Le traceur',iq,'appele ', &
+      &          trim(tnom_0(iq)),' est le fils de ',ipere,'appele ',trim(tnom_0(ipere))
+            nqfils(ipere)=nqfils(ipere)+1  
+            iqfils(nqfils(ipere),ipere)=iq
+            iqpere(iq)=ipere         
+            continu=.false.
+          else !if (tnom_transp(iq) == tnom_0(ipere)) then
+            ipere=ipere+1 
+            if (ipere.gt.nqtot) then
+                WRITE(lunout,*) 'Le traceur',iq,'appele ', &
+      &          trim(tnom_0(iq)),', est orpelin.'
+                CALL abort_gcm('infotrac_init','Un traceur est orphelin',1)
+            endif !if (ipere.gt.nqtot) then
+          endif !if (tnom_transp(iq) == tnom_0(ipere)) then
+        enddo !do while (continu)
+      endif !if (tnom_transp(iq) == 'air') then
+    enddo !DO iq=1,nqtot
+    WRITE(lunout,*) 'infotrac: nqperes=',nqperes    
+    WRITE(lunout,*) 'nqfils=',nqfils
+    WRITE(lunout,*) 'iqpere=',iqpere
+    WRITE(lunout,*) 'iqfils=',iqfils
+
+! Calculer le nombre de descendants à partir de iqfils et de nbfils
+    DO iq=1,nqtot    
+      generation=0
+      continu=.true.
+      ifils=iq
+      do while (continu)
+        ipere=iqpere(ifils)
+        if (ipere.gt.0) then
+         nqdesc(ipere)=nqdesc(ipere)+1   
+         nqdesc_tot=nqdesc_tot+1      
+         iqfils(nqdesc(ipere),ipere)=iq
+         ifils=ipere
+         generation=generation+1
+        else !if (ipere.gt.0) then
+         continu=.false.
+        endif !if (ipere.gt.0) then
+      enddo !do while (continu)    
+      WRITE(lunout,*) 'Le traceur ',iq,', appele ',trim(tnom_0(iq)),' est un traceur de generation: ',generation
+    enddo !DO iq=1,nqtot
+    WRITE(lunout,*) 'infotrac: nqdesc=',nqdesc
+    WRITE(lunout,*) 'iqfils=',iqfils
+    WRITE(lunout,*) 'nqdesc_tot=',nqdesc_tot
+
+! Interdire autres schémas que 10 pour les traceurs fils, et autres schémas
+! que 10 et 14 si des pères ont des fils
+    do iq=1,nqtot
+      if (iqpere(iq).gt.0) then
+        ! ce traceur a un père qui n'est pas l'air
+        ! Seul le schéma 10 est autorisé
+        if (iadv(iq)/=10) then
+           WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for sons'
+          CALL abort_gcm('infotrac_init','Sons should be advected by scheme 10',1)
+        endif
+        ! Le traceur père ne peut être advecté que par schéma 10 ou 14:
+        IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN
+          WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for fathers'
+          CALL abort_gcm('infotrac_init','Fathers should be advected by scheme 10 ou 14',1)
+        endif !IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN
+     endif !if (iqpere(iq).gt.0) the
+    enddo !do iq=1,nqtot
+
+
+! detecter quels sont les traceurs isotopiques parmi des traceurs
+    call infotrac_isoinit(tnom_0,nqtrue)
+        
 !-----------------------------------------------------------------------
 ! Finalize :
 !
-    DEALLOCATE(tnom_0, hadv, vadv)
+    DEALLOCATE(tnom_0, hadv, vadv,tnom_transp)
 
 
   END SUBROUTINE infotrac_init
 
+  SUBROUTINE infotrac_isoinit(tnom_0,nqtrue)
+
+#ifdef CPP_IOIPSL
+  use IOIPSL
+#else
+  ! if not using IOIPSL, we still need to use (a local version of) getin
+  use ioipsl_getincom
+#endif
+  implicit none
+ 
+    ! inputs
+    INTEGER nqtrue
+    CHARACTER(len=15) tnom_0(nqtrue)
+    
+    ! locals    
+    CHARACTER(len=3), DIMENSION(niso_possibles) :: tnom_iso
+    INTEGER, ALLOCATABLE,DIMENSION(:,:) :: nb_iso,nb_traciso
+    INTEGER, ALLOCATABLE,DIMENSION(:) :: nb_isoind
+    INTEGER :: ntraceurs_zone_prec,iq,phase,ixt,iiso,izone
+    CHARACTER(len=19) :: tnom_trac
+    INCLUDE "iniprint.h"
+
+    tnom_iso=(/'eau','HDO','O18','O17','HTO'/)
+
+    ALLOCATE(nb_iso(niso_possibles,nqo))
+    ALLOCATE(nb_isoind(nqo))
+    ALLOCATE(nb_traciso(niso_possibles,nqo))
+    ALLOCATE(iso_num(nqtot))
+    ALLOCATE(iso_indnum(nqtot))
+    ALLOCATE(zone_num(nqtot))
+    ALLOCATE(phase_num(nqtot))
+     
+    iso_num(:)=0
+    iso_indnum(:)=0
+    zone_num(:)=0
+    phase_num(:)=0
+    indnum_fn_num(:)=0
+    use_iso(:)=.false.  
+    nb_iso(:,:)=0  
+    nb_isoind(:)=0     
+    nb_traciso(:,:)=0
+    niso=0
+    ntraceurs_zone=0  
+    ntraceurs_zone_prec=0
+    ntraciso=0
+
+    do iq=nqo+1,nqtot
+       write(lunout,*) 'infotrac 569: iq,tnom_0(iq)=',iq,tnom_0(iq)
+       do phase=1,nqo    
+        do ixt= 1,niso_possibles   
+         tnom_trac=trim(tnom_0(phase))//'_'
+         tnom_trac=trim(tnom_trac)//trim(tnom_iso(ixt))
+         write(*,*) 'phase,ixt,tnom_trac=',phase,ixt,tnom_trac      
+         IF (tnom_0(iq) == tnom_trac) then
+          write(lunout,*) 'Ce traceur est un isotope'
+          nb_iso(ixt,phase)=nb_iso(ixt,phase)+1   
+          nb_isoind(phase)=nb_isoind(phase)+1   
+          iso_num(iq)=ixt
+          iso_indnum(iq)=nb_isoind(phase)
+          indnum_fn_num(ixt)=iso_indnum(iq)
+          phase_num(iq)=phase
+          write(lunout,*) 'iso_num(iq)=',iso_num(iq)
+          write(lunout,*) 'iso_indnum(iq)=',iso_indnum(iq)
+          write(lunout,*) 'indnum_fn_num(ixt)=',indnum_fn_num(ixt)
+          write(lunout,*) 'phase_num(iq)=',phase_num(iq)
+          goto 20
+         else if (iqpere(iq).gt.0) then          
+          if (tnom_0(iqpere(iq)) == tnom_trac) then
+           write(lunout,*) 'Ce traceur est le fils d''un isotope'
+           ! c'est un traceur d'isotope
+           nb_traciso(ixt,phase)=nb_traciso(ixt,phase)+1
+           iso_num(iq)=ixt
+           iso_indnum(iq)=indnum_fn_num(ixt)
+           zone_num(iq)=nb_traciso(ixt,phase)
+           phase_num(iq)=phase
+           write(lunout,*) 'iso_num(iq)=',iso_num(iq)
+           write(lunout,*) 'phase_num(iq)=',phase_num(iq)
+           write(lunout,*) 'zone_num(iq)=',zone_num(iq)
+           goto 20
+          endif !if (tnom_0(iqpere(iq)) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then
+         endif !IF (tnom_0(iq) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then
+        enddo !do ixt= niso_possibles
+       enddo !do phase=1,nqo
+  20   continue
+      enddo !do iq=1,nqtot
+
+      write(lunout,*) 'iso_num=',iso_num
+      write(lunout,*) 'iso_indnum=',iso_indnum
+      write(lunout,*) 'zone_num=',zone_num  
+      write(lunout,*) 'phase_num=',phase_num
+      write(lunout,*) 'indnum_fn_num=',indnum_fn_num
+
+      do ixt= 1,niso_possibles  
+
+        if (nb_iso(ixt,1).eq.1) then
+          ! on vérifie que toutes les phases ont le même nombre de
+          ! traceurs
+          do phase=2,nqo
+            if (nb_iso(ixt,phase).ne.nb_iso(ixt,1)) then
+              write(lunout,*) 'ixt,phase,nb_iso=',ixt,phase,nb_iso(ixt,phase)
+              CALL abort_gcm('infotrac_init','Phases must have same number of isotopes',1)
+            endif
+          enddo !do phase=2,nqo
+
+          niso=niso+1
+          use_iso(ixt)=.true.
+          ntraceurs_zone=nb_traciso(ixt,1)
+
+          ! on vérifie que toutes les phases ont le même nombre de
+          ! traceurs
+          do phase=2,nqo
+            if (nb_traciso(ixt,phase).ne.ntraceurs_zone) then
+              write(lunout,*) 'ixt,phase,nb_traciso=',ixt,phase,nb_traciso(ixt,phase)
+              write(lunout,*) 'ntraceurs_zone=',ntraceurs_zone
+              CALL abort_gcm('infotrac_init','Phases must have same number of tracers',1)
+            endif  
+          enddo  !do phase=2,nqo
+          ! on vérifie que tous les isotopes ont le même nombre de
+          ! traceurs
+          if (ntraceurs_zone_prec.gt.0) then               
+            if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
+              ntraceurs_zone_prec=ntraceurs_zone
+            else !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
+              write(*,*) 'ntraceurs_zone_prec,ntraceurs_zone=',ntraceurs_zone_prec,ntraceurs_zone   
+              CALL abort_gcm('infotrac_init', &
+               &'Isotope tracers are not well defined in traceur.def',1)           
+            endif !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
+           endif !if (ntraceurs_zone_prec.gt.0) then
+
+        else if (nb_iso(ixt,1).ne.0) then
+           WRITE(lunout,*) 'nqo,ixt=',nqo,ixt
+           WRITE(lunout,*) 'nb_iso(ixt,1)=',nb_iso(ixt,1)    
+           CALL abort_gcm('infotrac_init','Isotopes are not well defined in traceur.def',1)     
+        endif   !if (nb_iso(ixt,1).eq.1) then       
+    enddo ! do ixt= niso_possibles
+
+    ! dimensions isotopique:
+    ntraciso=niso*(ntraceurs_zone+1)
+    WRITE(lunout,*) 'niso=',niso
+    WRITE(lunout,*) 'ntraceurs_zone,ntraciso=',ntraceurs_zone,ntraciso    
+ 
+    ! flags isotopiques:
+    if (niso.gt.0) then
+        ok_isotopes=.true.
+    else
+        ok_isotopes=.false.
+    endif 
+    WRITE(lunout,*) 'ok_isotopes=',ok_isotopes
+ 
+    if (ok_isotopes) then
+        ok_iso_verif=.false.
+        call getin('ok_iso_verif',ok_iso_verif)
+        ok_init_iso=.false.
+        call getin('ok_init_iso',ok_init_iso)
+        tnat=(/1.0,155.76e-6,2005.2e-6,0.004/100.,0.0/)
+        alpha_ideal=(/1.0,1.01,1.006,1.003,1.0/)
+    endif !if (ok_isotopes) then  
+    WRITE(lunout,*) 'ok_iso_verif=',ok_iso_verif
+    WRITE(lunout,*) 'ok_init_iso=',ok_init_iso
+
+    if (ntraceurs_zone.gt.0) then
+        ok_isotrac=.true.
+    else
+        ok_isotrac=.false.
+    endif   
+    WRITE(lunout,*) 'ok_isotrac=',ok_isotrac
+
+    ! remplissage du tableau iqiso(ntraciso,phase)
+    ALLOCATE(iqiso(ntraciso,nqo))   
+    iqiso(:,:)=0     
+    do iq=1,nqtot
+        if (iso_num(iq).gt.0) then
+          ixt=iso_indnum(iq)+zone_num(iq)*niso
+          iqiso(ixt,phase_num(iq))=iq
+        endif
+    enddo
+    WRITE(lunout,*) 'iqiso=',iqiso
+
+    ! replissage du tableau index_trac(ntraceurs_zone,niso)
+    ALLOCATE(index_trac(ntraceurs_zone,niso))  
+    if (ok_isotrac) then
+        do iiso=1,niso
+          do izone=1,ntraceurs_zone 
+             index_trac(izone,iiso)=iiso+izone*niso
+          enddo
+        enddo
+    else !if (ok_isotrac) then      
+        index_trac(:,:)=0.0
+    endif !if (ok_isotrac) then
+    write(lunout,*) 'index_trac=',index_trac    
+
+! Finalize :
+    DEALLOCATE(nb_iso)
+
+  END SUBROUTINE infotrac_isoinit
+
 END MODULE infotrac