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

Timestamp:

Sep 22, 2021, 6:11:35 PM (3 years ago)

Author:

dcugnet

Message:

• fix of the delPhase function.
• getvar1 and getvar2 fixed and modified to avoid the usage of files with several time records and make the calls rather short.
• works again with iadv==0
• no more issues with tracers numbers (nqo, nqtot, etc.)
• fixes in the algebrical reduction routine used for "isotopes_parems.def" (containing simple expressions with variables that have to be substituted).
• still to be validated numerically

File:

• LMDZ6/branches/LMDZ-tracers/libf/dyn3dmem/dynetat0_loc.F90 (modified) (7 diffs)

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

@@ -7 +7 @@
 !-------------------------------------------------------------------------------
   USE parallel_lmdz
-  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 netcdf, ONLY: NF90_OPEN,  NF90_INQUIRE_DIMENSION, NF90_INQ_VARID,        &
+      NF90_NOWRITE, NF90_CLOSE, NF90_INQUIRE_VARIABLE,  NF90_GET_VAR,          &
+      NF90_GET_ATT, NF90_NoErr, NF90_INQUIRE
   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_ini, day_ref, itau_dyn, start_time
-  USE ener_mod,  ONLY: etot0, ptot0, ztot0, stot0, ang0
+  USE temps_mod, ONLY: annee_ref,day_ref,itau_dyn, &
+                       start_time,day_ini
+  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+  USE strings_mod, ONLY: strIdx
+  USE readTracFiles_mod, ONLY: known_phases, old_phases, nphases, phases_sep
 
   IMPLICIT NONE
@@ -41 +44 @@
   CHARACTER(LEN=256) :: sdum, var, modname, oldH2O
   INTEGER, PARAMETER :: length=100
-  INTEGER :: iq, fID, vID, idecal, ix, ip, ierr
+  INTEGER :: iq, fID, vID, idecal, ix, ip, ierr, ib, ie, nglo
   REAL    :: time, tab_cntrl(length)               !--- RUN PARAMS TABLE
   TYPE(tra), POINTER :: tr
@@ -122 +125 @@
   END IF
   CALL err(NF90_GET_VAR(fID,vID,time),"get",var)
-  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)
+  ib = ijb_v; ie = ije_v; nglo = ip1jm
+  CALL get_var2("vcov",     vcov(ib:ie,:), ib, ie, nglo)
+  ib = ijb_u; ie = ije_u; nglo = ip1jmp1
+  CALL get_var2("ucov",     ucov(ib:ie,:), ib, ie, nglo)
+  CALL get_var2("teta",     teta(ib:ie,:), ib, ie, nglo)
+  CALL get_var2("masse",   masse(ib:ie,:), ib, ie, nglo)
+  CALL get_var1("phisinit", phis(ib:ie),   ib, ie)
+  CALL get_var1("ps",         ps(ib:ie),   ib, ie)
 
 !--- Tracers
@@ -135 +140 @@
     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)
+      CALL get_var2(var, q(ib:ie,:,iq), ib, ie, nglo)
 #ifdef INCA
     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)
+      CALL get_var2('OX', q(ib:ie,:,iq), ib, ie, nglo)
 #endif
     ELSE IF(ix /= 0) THEN              !--- Old file, water: H2Ov/l/i instead of H2O-g/-l/-s
@@ -145 +150 @@
       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)
+        CALL get_var2(oldH2O, q(ib:ie,:,iq), ib, ie, nglo)
       END IF
     ELSE
       WRITE(lunout,*)TRIM(modname)//': Tracer <'//TRIM(var)//'> is missing => initialized to zero'
-      q(ijb_u:ije_u,:,iq)=0.
+      q(ib:ie,:,iq)=0.
       !--- CRisi: for isotopes, theoretical initialization using very simplified Rayleigh distillation law
       IF(niso > 0 .AND. tr%iso_num > 0) THEN
@@ -179 +184 @@
 
 
-SUBROUTINE get_var1(var, v, ib, ie, n_glo)
+SUBROUTINE get_var1(var, v, ib, ie)
+!--- Usable for fields up to rank 4 with single time record (last index)
+!--- Result: stacked in a vector. Used for 2D (single layer) fields.
   CHARACTER(LEN=*),  INTENT(IN)  :: var
   REAL,              INTENT(OUT) :: v(:)
-  INTEGER, OPTIONAL, INTENT(IN)  :: ib, ie, n_glo
+  INTEGER, OPTIONAL, INTENT(IN)  :: ib, ie
   REAL, ALLOCATABLE :: w(:,:,:,:), v_glo(:)
-  INTEGER :: nn(4), dids(4), k, nd, ntot
+  INTEGER :: n(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)))
+  n(:) = 1; DO k = 1, nd; ierr = NF90_INQUIRE_DIMENSION(fID, dids(k), len=n(k)); END DO
+  IF(is_rec(fID, dids(nd)) .AND. n(nd) /= 1) &
+    CALL abort_gcm(TRIM(modname), 'Several records records for <'//TRIM(var)//'>')
+  ntot = PRODUCT(n(1:nd))
+  ALLOCATE(w(n(1), n(2), n(3), n(4)), v_glo(ntot))
   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
+  v_glo(:) = RESHAPE(w, [ntot]); DEALLOCATE(w)
+  IF(PRESENT(ib).AND.PRESENT(ie)) THEN; v(:) = v_glo(ib:ie); ELSE; v(:) = v_glo(:); END IF
   DEALLOCATE(v_glo)
 END SUBROUTINE get_var1
@@ -203 +207 @@
 
 SUBROUTINE get_var2(var, v, ib, ie, n_glo)
+!--- Usable for fields up to rank 4 with one or several time records (last index)
+!--- Result: stacked in a 2D array (1st/2nd index: horizontal/vertical). Used for 3D (several layers) fields.
   CHARACTER(LEN=*), INTENT(IN)  :: var
   REAL,             INTENT(OUT) :: v(:,:)
   INTEGER,          INTENT(IN)  :: ib, ie, n_glo
   REAL, ALLOCATABLE :: w(:,:,:,:), v_glo(:,:)
-  INTEGER :: nn(4), dids(4), k, nd, nh, nv, tid
+  INTEGER :: n(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
+  n(:) = 1; DO k = 1, nd; ierr = NF90_INQUIRE_DIMENSION(fID, dids(k), len=n(k)); END DO
+  IF(is_rec(fID, dids(nd))) THEN
+    IF(n(nd) /= 1)  CALL abort_gcm(TRIM(modname), 'Several records records for <'//TRIM(var)//'>.')
+    nh = PRODUCT(n(1:nd-2)); nv = n(nd-1)
   ELSE
-    nh = PRODUCT(nn(1:nd-1)); nv = nn(nd)
-  END IF
-  ALLOCATE(w(nn(1), nn(2), nn(3), nn(4)))
+    nh = PRODUCT(n(1:nd  )); nv = n(nd)
+  END IF
+  IF(nh/=n_glo .OR. nv/=llm) CALL abort_gcm(TRIM(modname), 'Shape mismatch for "'//TRIM(var)//'"')
+  ALLOCATE(w(n(1), n(2), n(3), n(4)), v_glo(nh,nv))
   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,:)
+  v_glo(:,:) = RESHAPE(w, [nh, nv]); DEALLOCATE(w)
+  v(:,:) = v_glo(ib:ie,:)
   DEALLOCATE(v_glo)
 END SUBROUTINE get_var2
+
+
+LOGICAL FUNCTION is_rec(fID, did) RESULT(lrec)
+!--- Check whether the file has a record dimension, detected as UNLIMITED diemnsion or using the attribute "units".
+  INTEGER, INTENT(IN) :: fID, did
+  INTEGER :: vid
+  CHARACTER(LEN=256) :: recn, ratt
+  !--- Check the "units" attribute of the last dimensional variable to detect record axis.
+  lrec = NF90_INQUIRE_DIMENSION  (fID, did, name=recn)    == NF90_NOERR
+  IF(lrec) lrec = NF90_INQ_VARID (fID, recn, vid)         == NF90_NOERR
+  IF(lrec) lrec = NF90_GET_ATT   (fID, vid, "units", ratt)== NF90_NOERR
+  IF(lrec) lrec = INDEX(ratt, " since ") /= 0
+END FUNCTION is_rec