Changeset 3957 for LMDZ6/branches/LMDZ-tracers/libf/dyn3dmem

Timestamp:

Jul 11, 2021, 11:39:01 PM (3 years ago)

Author:

dcugnet

Message:

In readTracFiles: the separator between the tracer name and its phase is no longer hardcoded and equal to "-",
but is a parameter ("phases_sep") which default value is "_".
Few more fixes.

Location:

LMDZ6/branches/LMDZ-tracers/libf/dyn3dmem

Files:

• dynetat0_loc.F90 (modified) (6 diffs)
• vlspltqs_loc.F (modified) (1 diff)

LMDZ6/branches/LMDZ-tracers/libf/dyn3dmem/dynetat0_loc.F90

@@ -7 +7 @@
 !-------------------------------------------------------------------------------
   USE parallel_lmdz
-  USE infotrac, ONLY: nqtot, niso, tracers, iTraPha, tnat, alpha_ideal, tra
-  USE netcdf, ONLY: NF90_OPEN,  NF90_INQUIRE_DIMENSION, NF90_INQ_VARID,        &
-      NF90_NOWRITE, NF90_CLOSE, NF90_INQUIRE_VARIABLE,  NF90_GET_VAR, NF90_NoErr
+  USE readTracFiles_mod, ONLY: known_phases, old_phases, nphases, phases_sep
+  USE infotrac,    ONLY: nqtot, niso, tracers, iTraPha, tnat, alpha_ideal, tra
+  USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQ_VARID, NF90_INQUIRE_DIMENSION, &
+           NF90_INQUIRE, NF90_CLOSE, NF90_GET_VAR, NF90_NoErr,     NF90_INQUIRE_VARIABLE
+  USE strings_mod, ONLY: strIdx
   USE control_mod, ONLY: planet_type
   USE assert_eq_m, ONLY: assert_eq
   USE comvert_mod, ONLY: pa,preff
-  USE comconst_mod, ONLY: cpp, daysec, dtvr, g, im, jm, kappa, lllm, &
-                          omeg, rad
+  USE comconst_mod, ONLY: cpp, daysec, dtvr, g, im, jm, kappa, lllm, omeg, rad
   USE logic_mod, ONLY: fxyhypb, ysinus
   USE serre_mod, ONLY: clon, clat, grossismx, grossismy
-  USE temps_mod, ONLY: annee_ref,day_ref,itau_dyn, &
-                       start_time,day_ini,hour_ini
-  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
-  
+  USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn, start_time
+  USE ener_mod,  ONLY: etot0, ptot0, ztot0, stot0, ang0
+
   IMPLICIT NONE
   include "dimensions.h"
@@ -39 +39 @@
 !===============================================================================
 ! Local variables:
-  CHARACTER(LEN=256) :: sdum, var, modname
+  CHARACTER(LEN=256) :: sdum, var, modname, oldH2O
   INTEGER, PARAMETER :: length=100
-  INTEGER :: iq, fID, vID, idecal, ierr
+  INTEGER :: iq, fID, vID, idecal, ix, ip, ierr
   REAL    :: time, tab_cntrl(length)               !--- RUN PARAMS TABLE
-  REAL,             ALLOCATABLE :: vcov_glo(:,:),masse_glo(:,:),   ps_glo(:)
-  REAL,             ALLOCATABLE :: ucov_glo(:,:),    q_glo(:,:), phis_glo(:)
-  REAL,             ALLOCATABLE :: teta_glo(:,:)
   TYPE(tra), POINTER :: tr
 !-------------------------------------------------------------------------------
@@ -110 +107 @@
   CALL check_dim(jm,jjm,'jm','jm')
   CALL check_dim(lllm,llm,'lm','lllm')
-  CALL get_var1("rlonu",rlonu)
-  CALL get_var1("rlatu",rlatu)
-  CALL get_var1("rlonv",rlonv)
-  CALL get_var1("rlatv",rlatv)
-  CALL get_var1("cu"  ,cu)
-  CALL get_var1("cv"  ,cv)
-  CALL get_var1("aire",aire)
+  CALL get_var1("rlonu", rlonu)
+  CALL get_var1("rlatu", rlatu)
+  CALL get_var1("rlonv", rlonv)
+  CALL get_var1("rlatv", rlatv)
+  CALL get_var1("cu",       cu)
+  CALL get_var1("cv",       cv)
+  CALL get_var1("aire",   aire)
 
   var="temps"
@@ -125 +122 @@
   END IF
   CALL err(NF90_GET_VAR(fID,vID,time),"get",var)
-
-  ALLOCATE(phis_glo(ip1jmp1))
-  CALL get_var1("phisinit",phis_glo)
-  phis (ijb_u:ije_u)  =phis_glo(ijb_u:ije_u);    DEALLOCATE(phis_glo)
-
-  ALLOCATE(ucov_glo(ip1jmp1,llm))
-  CALL get_var2("ucov",ucov_glo)
-  ucov (ijb_u:ije_u,:)=ucov_glo(ijb_u:ije_u,:);  DEALLOCATE(ucov_glo)
-
-  ALLOCATE(vcov_glo(ip1jm,llm))
-  CALL get_var2("vcov",vcov_glo)
-  vcov (ijb_v:ije_v,:)=vcov_glo(ijb_v:ije_v,:);  DEALLOCATE(vcov_glo)
-
-  ALLOCATE(teta_glo(ip1jmp1,llm))
-  CALL get_var2("teta",teta_glo)
-  teta (ijb_u:ije_u,:)=teta_glo(ijb_u:ije_u,:);  DEALLOCATE(teta_glo)
-
-  ALLOCATE(masse_glo(ip1jmp1,llm))
-  CALL get_var2("masse",masse_glo)
-  masse(ijb_u:ije_u,:)=masse_glo(ijb_u:ije_u,:); DEALLOCATE(masse_glo)
-  
-  ALLOCATE(ps_glo(ip1jmp1))
-  CALL get_var1("ps",ps_glo)
-  ps   (ijb_u:ije_u)  =   ps_glo(ijb_u:ije_u);   DEALLOCATE(ps_glo)
+  CALL get_var1("phisinit", phis, ijb_u, ije_u, ip1jmp1)
+  CALL get_var2("ucov",     ucov, ijb_u, ije_u, ip1jmp1)
+  CALL get_var2("vcov",     vcov, ijb_v, ije_v, ip1jm)
+  CALL get_var2("teta",     teta, ijb_u, ije_u, ip1jmp1)
+  CALL get_var2("masse",   masse, ijb_u, ije_u, ip1jmp1)
+  CALL get_var1("ps",         ps, ijb_u, ije_u, ip1jmp1)
 
 !--- Tracers
-  ALLOCATE(q_glo(ip1jmp1,llm))
   DO iq=1,nqtot
     tr => tracers(iq)
     var = tr%name
-    IF(NF90_INQ_VARID(fID,var,vID)==NF90_NoErr) THEN
-      CALL get_var2(var ,q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:); CYCLE
+    ix = strIdx([('H2O'//phases_sep//known_phases(ip:ip), ip=1, nphases)], var)
+    IF(NF90_INQ_VARID(fID, var, vID) == NF90_NoErr) THEN
+      CALL get_var2(var, q(:,:,iq), ijb_u, ije_u, ip1jmp1)
 #ifdef INCA
-    ELSE IF(NF90_INQ_VARID(fID,"OX",vID) == NF90_NoErr .AND. var == 'O3') THEN
-      WRITE(lunout,*) 'Tracer O3 is missing - it is initialized to OX'
-      CALL get_var2("OX",q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:); CYCLE
+    ELSE IF(NF90_INQ_VARID(fID, 'OX', vID) == NF90_NoErr .AND. var == 'O3') THEN
+      WRITE(lunout,*)TRIM(modname)//': Tracer <O3> is missing => initialized to OX'
+      CALL get_var2('OX', q(:,:,iq), ijb_u, ije_u, ip1jmp1)
 #endif
-    END IF
-    WRITE(lunout,*)TRIM(modname)//": Tracer <"//TRIM(var)//"> is missing"
-    WRITE(lunout,*)"         It is hence initialized to zero"
-    q(ijb_u:ije_u,:,iq)=0.
-   !--- CRisi: for isotops, theoretical initialization using very simplified
-   !           Rayleigh distillation las.
-    IF(niso > 0 .AND. tr%iso_num > 0) THEN
-      IF(tr%iso_zon == 0) q(:,:,iq) = q(:,:,tr%iprnt)         *        tnat(tr%iso_num) &
-                                   * (q(:,:,tr%iprnt)/30.e-3)**(alpha_ideal(tr%iso_num)-1)
-      IF(tr%iso_zon == 1) q(:,:,iq) = q(:,:,iTraPha(tr%iso_num,tr%iso_pha))
+    ELSE IF(ix /= 0) THEN              !--- Old file, water: H2Ov/l/i instead of H2O-g/-l/-s
+      oldH2O = 'H2O'//old_phases(ix:ix)
+      IF(NF90_INQ_VARID(fID, oldH2O, vID) == NF90_NoErr) THEN
+        WRITE(lunout,*)TRIM(modname)//': Tracer <'//TRIM(var)//'> is missing => initialized to '//TRIM(oldH2O)
+        CALL get_var2(oldH2O, q(:,:,iq), ijb_u, ije_u, ip1jmp1)
+      END IF
+    ELSE
+      WRITE(lunout,*)TRIM(modname)//': Tracer <'//TRIM(var)//'> is missing => initialized to zero'
+      q(ijb_u:ije_u,:,iq)=0.
+      !--- CRisi: for isotopes, theoretical initialization using very simplified Rayleigh distillation law
+      IF(niso > 0 .AND. tr%iso_num > 0) THEN
+        IF(tr%iso_zon == 0) q(:,:,iq) = q(:,:,tr%iprnt)         *        tnat(tr%iso_num) &
+                                      *(q(:,:,tr%iprnt)/30.e-3)**(alpha_ideal(tr%iso_num)-1)
+        IF(tr%iso_zon == 1) q(:,:,iq) = q(:,:,iTraPha(tr%iso_num,tr%iso_pha))
+      END IF
     END IF
   END DO
-  DEALLOCATE(q_glo)
   CALL err(NF90_CLOSE(fID),"close",fichnom)
   day_ini=day_ini+INT(time)
@@ -190 +173 @@
     s1='value of '//TRIM(str1)//' ='
     s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
-    WRITE(sdum,'(10x,a,i4,2x,a,i4)'),s1,n1,s2,n2
+    WRITE(sdum,'(10x,a,i4,2x,a,i4)'),TRIM(s1),n1,TRIM(s2),n2
     CALL ABORT_gcm(TRIM(modname),TRIM(sdum),1)
   END IF
@@ -196 +179 @@
 
 
-SUBROUTINE get_var1(var,v)
-  CHARACTER(LEN=*), INTENT(IN)  :: var
-  REAL,             INTENT(OUT) :: v(:)
-  REAL,             ALLOCATABLE :: w2(:,:), w3(:,:,:)
-  INTEGER :: nn(3), dids(3), k, nd, ntot
-
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  ierr=NF90_INQUIRE_VARIABLE(fID,vID,ndims=nd)
-  IF(nd==1) THEN
-    CALL err(NF90_GET_VAR(fID,vID,v),"get",var); RETURN
-  END IF
-  ierr=NF90_INQUIRE_VARIABLE(fID,vID,dimids=dids)
-  DO k=1,nd; ierr=NF90_INQUIRE_DIMENSION(fID,dids(k),len=nn(k)); END DO
-  ntot=PRODUCT(nn(1:nd))
-  SELECT CASE(nd)
-    CASE(2); ALLOCATE(w2(nn(1),nn(2)))
-      CALL err(NF90_GET_VAR(fID,vID,w2),"get",var)
-      v=RESHAPE(w2,[ntot]); DEALLOCATE(w2)
-    CASE(3); ALLOCATE(w3(nn(1),nn(2),nn(3)))
-      CALL err(NF90_GET_VAR(fID,vID,w3),"get",var)
-      v=RESHAPE(w3,[ntot]); DEALLOCATE(w3)
-  END SELECT
+SUBROUTINE get_var1(var, v, ib, ie, n_glo)
+  CHARACTER(LEN=*),  INTENT(IN)  :: var
+  REAL,              INTENT(OUT) :: v(:)
+  INTEGER, OPTIONAL, INTENT(IN)  :: ib, ie, n_glo
+  REAL, ALLOCATABLE :: w(:,:,:,:), v_glo(:)
+  INTEGER :: nn(4), dids(4), k, nd, ntot
+  CALL err(NF90_INQ_VARID(fID, var, vID), "inq", var)
+  ierr = NF90_INQUIRE_VARIABLE(fID, vID, dimids=dids, ndims=nd)
+  nn(:) = 1; DO k=1,nd; ierr = NF90_INQUIRE_DIMENSION(fID, dids(k), len=nn(k)); END DO
+  ntot = PRODUCT(nn(1:nd))
+  ALLOCATE(w(nn(1), nn(2), nn(3), nn(4)))
+  CALL err(NF90_GET_VAR(fID, vID, w), "get", var)
+  ALLOCATE(v_glo(ntot)); v_glo=RESHAPE(w, [ntot]); DEALLOCATE(w)
+  IF(PRESENT(n_glo).AND.PRESENT(ib).AND.PRESENT(ie)) THEN
+    IF(ntot/=n_glo) CALL abort_gcm(TRIM(modname), 'Shape mismatch for "'//TRIM(var)//'"')
+    v(ib:ie) = v_glo(ib:ie)
+  ELSE
+    v(:) = v_glo(:)
+  END IF
+  DEALLOCATE(v_glo)
 END SUBROUTINE get_var1
 
 
-SUBROUTINE get_var2(var,v)
+SUBROUTINE get_var2(var, v, ib, ie, n_glo)
   CHARACTER(LEN=*), INTENT(IN)  :: var
   REAL,             INTENT(OUT) :: v(:,:)
-  REAL,             ALLOCATABLE :: w4(:,:,:,:)
-  INTEGER :: nn(4), dids(4), k, nd
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  ierr=NF90_INQUIRE_VARIABLE(fID,vID,dimids=dids,ndims=nd)
-  DO k=1,nd; ierr=NF90_INQUIRE_DIMENSION(fID,dids(k),len=nn(k)); END DO
-  ALLOCATE(w4(nn(1),nn(2),nn(3),nn(4)))
-  CALL err(NF90_GET_VAR(fID,vID,w4),"get",var)
-  v=RESHAPE(w4,[nn(1)*nn(2),nn(3)]); DEALLOCATE(w4)
+  INTEGER,          INTENT(IN)  :: ib, ie, n_glo
+  REAL, ALLOCATABLE :: w(:,:,:,:), v_glo(:,:)
+  INTEGER :: nn(4), dids(4), k, nd, nh, nv, tid
+  CALL err(NF90_INQ_VARID(fID, var, vID), "inq", var)
+  ierr = NF90_INQUIRE_VARIABLE(fID, vID, dimids=dids, ndims=nd)
+  nn(:) = 1; DO k=1,nd; ierr = NF90_INQUIRE_DIMENSION(fID, dids(k), len=nn(k)); END DO
+  IF(NF90_INQUIRE(fID, unlimitedDimId=tid) == NF90_NOERR) THEN
+    nh = PRODUCT(nn(1:nd-2)); nv = nn(nd-1); nn(nd) = 1
+  ELSE
+    nh = PRODUCT(nn(1:nd-1)); nv = nn(nd)
+  END IF
+  ALLOCATE(w(nn(1), nn(2), nn(3), nn(4)))
+  CALL err(NF90_GET_VAR(fID, vID, w), "get", var)
+  ALLOCATE(v_glo(nh, nv)); v_glo = RESHAPE(w, [nh, nv]); DEALLOCATE(w)
+  IF(nh/=n_glo .OR. nv/=llm) CALL abort_gcm(TRIM(modname), 'Shape mismatch for "'//TRIM(var)//'"')
+  v(ib:ie,:) = v_glo(ib:ie,:)
+  DEALLOCATE(v_glo)
 END SUBROUTINE get_var2
 
 
-SUBROUTINE err(ierr,typ,nam)
+SUBROUTINE err(ierr, typ, nam)
   INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
   CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION

LMDZ6/branches/LMDZ-tracers/libf/dyn3dmem/vlspltqs_loc.F

@@ -481 +481 @@
       INTEGER ijbm,ijem
 
+      tr => tracers(iq)
       ijb=ij_begin-2*iip1
       ije=ij_end+2*iip1