Changeset 136 for trunk/SOURCES

Timestamp:

10/13/17 18:21:01 (7 years ago)

Author:

dumas

Message:

Flux output and short netcdf output for initMIP

Location:

trunk/SOURCES

Files:

• 3D-physique-gen_mod.f90 (modified) (1 diff)
• Ant40_files/output_anta40_mod-0.4.f90 (modified) (8 diffs)
• Makefile.grisli.inc (modified) (2 diffs)
• ablation_bord.f90 (modified) (1 diff)
• bilan_eau_mod.f90 (modified) (3 diffs)
• initial-0.3.f90 (modified) (2 diffs)
• steps_time_loop.f90 (modified) (3 diffs)

trunk/SOURCES/3D-physique-gen_mod.f90

@@ -226 +226 @@
   integer,dimension(nx,ny) :: FRONTFACEY!< type de front sur les faces y
   integer,dimension(nx,ny) :: gr_line_schoof ! points ou on impose le flux de schoof (pour sorties)
+  integer,dimension(nx,ny) :: gr_line   !< points grounding line pour les sorties
   integer,dimension(nx,ny) :: ICE       !< presence de glace au point considere,
   !< seuil 0 si pose, 1 si flottant

trunk/SOURCES/Ant40_files/output_anta40_mod-0.4.f90

@@ -15 +15 @@
 USE module3D_phy
 use bilan_eau_mod
-
+use netcdf
+use io_netcdf_grisli
+use bilan_flux_mod
 
 implicit none
       
 !real ::  vol ; !integer :: np
+integer :: nflot                      !< nbr de point flottant
 real ::  bmean                        !<
 real ::  accmean                      !< accumulation moyenne
@@ -35 +38 @@
 real ::  hdotmean                     !< moyenne derivee / temps de h
 real ::  bdotmean                     !< moyenne bedrock derive / temps
-real :: volf                          !< volume au dessus de la flottaison
+real ::  volf                         !< volume au dessus de la flottaison
+
+real :: lim                           !< Total ice mass
+real :: iareag                        !< surface posee
+real :: iareaf                        !< surface flottante
+real :: tendacabf                     !< Total SMB flux
+real :: tendlibmassbf                 !< Total Basal mass balance flux
+real :: tendlicalvf                   !< Total calving flux
+real :: tendligroundf                 !< Total grounding line flux
+
+real,dimension(nx,ny) :: corrsurf     !< facteur de correction de la surface
+real, parameter :: ice_density=910   !< densite de la glace pour conversion en masse
+
+! variables netcdf
+integer,parameter :: ncshortout=1 ! 1 sorties netcdf short initMIP
+integer,parameter :: nvar=8 ! nombre de variables dans le fichier de sortie temporel netcdf 
+integer :: ncid
+integer :: status
+integer :: timeDimID
+integer :: timeVarID
+integer,dimension(nvar) :: varID
+integer :: nbtimeout  ! index time output
+
+real,dimension(nvar) :: var_shortoutput
 
 
@@ -41 +67 @@
 
 subroutine init_outshort
+
+double precision,dimension(:,:),pointer      :: tab               !< tableau 2d real pointer
+character(len=100),dimension(nvar) :: namevar ! name, standard_name, long_name, units
+character(len=100),dimension(nvar) :: standard_name
+character(len=100),dimension(nvar) :: long_name
+character(len=100),dimension(nvar) :: units
 
 !ndisp sorite courte tous les ndisp
 NDISP=100
+
+
+
+if (ncshortout.eq.1) then  ! ecriture netcdf
+
+! lecture du fichier avec les corrections de surface
+  call Read_Ncdf_var('z',trim(DIRNAMEINP)//'/corrsurf-initMIP-16km.grd',tab)
+  corrsurf(:,:)=tab(:,:)
+
+  open(568,file=trim(dirsource)//'/Fichiers-parametres/short-initMIPnc.dat',status='old')
+! lecture en-tete
+  read(568,*)
+  read(568,*)
+  read(568,*)
+  read(568,*)
+  read(568,*)
+! lecture des infos sur les variables :
+  do k=1,nvar
+    read(568,'(a100)') namevar(k)
+    print*,trim(namevar(k))
+    read(568,'(a100)') standard_name(k)
+    print*,trim(standard_name(k))
+    read(568,'(a100)') long_name(k)
+    print*,trim(long_name(k))
+    read(568,'(a100)') units(k)
+    print*,trim(units(k))
+    read(568,*)
+  enddo
+  close(568)
+  print*,'apres lecture fichier descriptions variables'
+! Fichier Netcdf initMIP
+! creation du fichier Netcdf :
+  status=nf90_create(path = 'short'//runname//'.nc', cmode = nf90_clobber, ncid = ncid)
+  if (status /= nf90_noerr) call handle_err(status)
+  print*,'apres nf90_create'
+
+! definition des dimension :
+  status=nf90_def_dim(ncid, name="time", len=NF90_UNLIMITED, dimid=timeDimID)
+  if (status /= nf90_noerr) call handle_err(status)
+  print*,'apres 2'
+  status=nf90_def_var(ncid, name="time", xtype=nf90_float, dimids=(/ timeDimID/), varid=timeVarID) 
+  if (status /= nf90_noerr) call handle_err(status)
+  status=nf90_put_att(ncid, timeVarID, "standard_name", "time")
+  if (status /= nf90_noerr) call handle_err(status)
+  print*,'apres 3'
+  status=nf90_put_att(ncid, timeVarID,"units", "years since 2007-01-01 00:00:00")
+  if (status /= nf90_noerr) call handle_err(status)
+  print*,'apres 4'
+
+! definition des variables de sortie :
+  do k=1,nvar  ! boucle sur le nbr de variable a definir
+    status=nf90_def_var(ncid, name=trim(namevar(k)), xtype=nf90_float, dimids= &
+          (/ timeDimID /), varid=varID(k))
+  print*,'apres 5'
+    if (status /= nf90_noerr) call handle_err(status)
+    status=nf90_put_att(ncid, varID(k), "standard_name", trim(standard_name(k)))
+    print*,'apres 6',trim(standard_name(k))
+    if (status /= nf90_noerr) call handle_err(status)
+    status=nf90_put_att(ncid, varID(k), "long_name", trim(long_name(k)))
+    print*,'apres 7',trim(long_name(k))
+    if (status /= nf90_noerr) call handle_err(status)
+    status=nf90_put_att(ncid, varID(k), "units", trim(units(k)))
+    print*,'apres 8', trim(units(k))
+    if (status /= nf90_noerr) call handle_err(status)
+  enddo
+  print*,'apres boucle'
+
+! fin de la definition du fchier :
+  status=nf90_enddef(ncid)
+  if (status /= nf90_noerr) call handle_err(status)
+  nbtimeout = 0 ! initialisation compteur sorties axe time
+else ! pas de sortie netcdf et sans correction de surface
+  corrsurf(:,:)=1.
+endif
+
 end subroutine init_outshort
 
@@ -54 +161 @@
 !------------------
 real ::  smax
+
       vol=0. 
       np=0
+      nflot=0
       hmax=0. 
       smax=0.
@@ -74 +183 @@
       bdotmean=0.
       volf=0.
- 
+      lim=0.
+      tendacabf=0.
+      iareag=0.
+      iareaf=0.
+      tendlicalvf=0.
+      tendligroundf=0.
+      tendlibmassbf=0.
 ! 2_preparing outputs
 !--------------------      
-    do i=1,nx 
-      do j=1,ny
-        if (.not.flot(i,j)) then 
-!       if (h(i,j).gt.1.) then 
-          np=np+1
-          vol=vol+h(i,j) 
-
-!        calcul de la hauteur au dessus de la flottaison
-         if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
-               volf=volf+h(i,j)
-         else
-            volf=volf+h(i,j)-row/ro*(sealevel-b(i,j))
-         endif
-
-
-          if (h(i,j).gt.hmax) hmax=h(i,j) 
-          if (s(i,j).gt.smax) smax=s(i,j) 
-          bmean=bm(i,j)+bmean 
-          accmean=acc(i,j)+accmean
-          tbmean=tbmean+t(i,j,nz)
-          tbdotmean=tbdotmean+tbdot(i,j)
-          vsmean=vsmean+sqrt(ux(i,j,1)**2+uy(i,j,1)**2)
+    do j=1,ny
+      do i=1,nx
+        if (ice(i,j).eq.1) then ! point englace
+          if (.not.flot(i,j)) then ! point pose 
+            np=np+1
+            vol=vol+h(i,j)
+            iareag=iareag+1.*corrsurf(i,j)                            ! surface englacee posee
+
+!         calcul de la hauteur au dessus de la flottaison
+            if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
+              volf=volf+h(i,j)*corrsurf(i,j)                          ! volume au-dessus de la flottaison
+            else
+              volf=volf+(h(i,j)-row/ro*(sealevel-b(i,j)))*corrsurf(i,j) ! volume au-dessus de la flottaison
+            endif
+
+            if (h(i,j).gt.hmax) hmax=h(i,j) 
+            if (s(i,j).gt.smax) smax=s(i,j) 
+            bmean=bm(i,j)+bmean 
+            accmean=acc(i,j)+accmean
+            tbmean=tbmean+t(i,j,nz)
+            tbdotmean=tbdotmean+tbdot(i,j)
+            vsmean=vsmean+sqrt(ux(i,j,1)**2+uy(i,j,1)**2)
 !          vsdotmean=vsdotmean+vsdot(i,j)
 !          uzsmean=uzsmean+uz(i,j,1)
-          uzsdotmean=uzsdotmean+uzsdot(i,j)
-          uzkmean=uzkmean+uzk(i,j)
-          hdotmean=hdotmean+abs(hdot(i,j)) 
-          bdotmean=bdotmean+abs(bdot(i,j)) 
-          bmeltmean=bmeltmean+bmelt(i,j)
+            uzsdotmean=uzsdotmean+uzsdot(i,j)
+            uzkmean=uzkmean+uzk(i,j)
+            hdotmean=hdotmean+abs(hdot(i,j)) 
+            bdotmean=bdotmean+abs(bdot(i,j)) 
+            bmeltmean=bmeltmean+bmelt(i,j)
+          else ! point flottant
+            iareaf=iareaf+1.*corrsurf(i,j)                               ! surface flottante
+          endif
+          lim=lim+h(i,j)*corrsurf(i,j)                                   ! volume total de glace
+          tendacabf=tendacabf+bm_dtt(i,j)*corrsurf(i,j)/dtt              ! smb surface
+          tendlibmassbf=tendlibmassbf-bmelt_dtt(i,j)*corrsurf(i,j)/dtt   ! fonte basale     
         endif
-        calvmean=calvmean+calv(i,j)  
-        ablbordmean=ablbordmean+ablbord(i,j)
+        tendlicalvf=tendlicalvf+calv_dtt(i,j)*corrsurf(i,j)/dtt        ! calving
+        tendlibmassbf=tendlibmassbf-ablbord_dtt(i,j)*corrsurf(i,j)/dtt ! partie ablbord de la fonte basale
+        tendligroundf=tendligroundf+grline_dtt(i,j)*corrsurf(i,j)/dtt  ! flux a la grounding line
       end do
-      end do
+    end do
 
 
@@ -116 +238 @@
         hmean=vol/np 
         vol=vol*dx*dy 
-        volf=volf*dx*dy
+        volf=volf*dx*dy*ice_density
         bmean=bmean/np 
         accmean=accmean/np 
@@ -130 +252 @@
         uzsdotmean=uzsdotmean/np
         uzkmean=uzkmean/np
-        hdotmean=hdotmean/np 
-      endif 
-
+        hdotmean=hdotmean/np
+      endif
+      lim=lim*dx*dy*ice_density
+      iareag=iareag*dx*dy 
+      iareaf=iareaf*dx*dy
+      tendacabf=tendacabf*dx*dy*ice_density/86400.
+      tendlibmassbf=tendlibmassbf*ice_density/86400.
+      tendlicalvf=tendlicalvf*ice_density/86400.
+      tendligroundf=tendligroundf*ice_density/86400.
+      
+      
       bdotmean=bdotmean/nx/ny  
 
@@ -154 +284 @@
 !940   format('%%%% ',a,'   time=',f8.0,' %%%%')
 
+print*, 'dans shortoutput avant netcdf ncshortout', ncshortout
+
+if (ncshortout.eq.1) then  ! ecriture netcdf
+  ! Total ice mass
+  var_shortoutput(1)=lim
+  ! Mass above floatation
+  var_shortoutput(2)=volf
+  ! Grounded ice area
+  var_shortoutput(3)=iareag
+  ! Floating ice area
+  var_shortoutput(4)=iareaf
+  ! Total SMB flux 
+  var_shortoutput(5)=tendacabf
+  ! Total Basal mass balance flux
+  var_shortoutput(6)=tendlibmassbf
+  ! Total calving flux 
+  var_shortoutput(7)=tendlicalvf
+  ! Total grounding line flux
+  var_shortoutput(8)=tendligroundf
+
+  nbtimeout=nbtimeout+1
+  print*,'ecriture shortnc', nbtimeout
+  status=nf90_put_var(ncid, timeVarID, time, start=(/nbtimeout/))
+  if (status /= nf90_noerr) call handle_err(status)
+
+  do k=1,nvar  ! boucle sur le nbr de variable a ecrire
+    status=nf90_put_var(ncid, varID(k), var_shortoutput(k),start=(/nbtimeout/))
+    if (status /= nf90_noerr) call handle_err(status)
+  enddo
+  status=nf90_sync(ncid)
+  if (status /= nf90_noerr) call handle_err(status)
+endif
+
+
 end subroutine shortoutput
+
+subroutine handle_err(status)
+  integer, intent(in) :: status
+  if (status /= nf90_noerr) then
+     print*,trim(nf90_strerror(status))
+     stop "stopped"
+  end if
+end subroutine handle_err
+
 end module  output_antarcti_mod

trunk/SOURCES/Makefile.grisli.inc

@@ -77 +77 @@
         calving_frange.o no_calving.o no_lakes.o \
         out_profile_mod.o printtable_mod.o mix-SIA-L1_mod.o \
-        furst_schoof_mod.o \
+        furst_schoof_mod.o bilan_flux_output_mod.o \
         relaxation_water_diffusion.o \
         prescribe-H-i2s_mod.o  \
@@ -97 +97 @@
         calving_frange.o no_calving.o no_lakes.o \
         out_profile_mod.o printtable_mod.o mix-SIA-L1_mod.o \
-        furst_schoof_mod.o \
+        furst_schoof_mod.o bilan_flux_output_mod.o \
         relaxation_water_diffusion.o \
         prescribe-H-i2s_mod.o   \

trunk/SOURCES/ablation_bord.f90

@@ -55 +55 @@
 ablbord_dtt(:,:) = ablbord_dtt(:,:) + ablbord(:,:)
 
-ablbord(:,:)=0.
+!cdc deplace dans bilan_flux_output_mod ablbord(:,:)=0.
 
 end subroutine ablation_bord

trunk/SOURCES/bilan_eau_mod.f90

@@ -31 +31 @@
 real,dimension(nx,ny) :: Bmelt_dtt    !< basal melting on ice points accumulated during dtt
 real,dimension(nx,ny) :: calv_dtt     !< calving sur dtt (pour calcul bilan d'eau)
+real,dimension(nx,ny) :: archimtab    ! point pose si > 0
+real,dimension(nx,ny) :: grline_dtt   ! grounding line flux during dtt
 
 real :: sum_H_old
 real :: diff_H
+real :: alpha_flot                      ! ro/row 
 
 real,dimension(nx,ny) :: bm_dt,bmelt_dt
@@ -47 +50 @@
         bm_dt(:,:)=0.
         bmelt_dt(:,:)=0.
+  alpha_flot=ro/row
 end subroutine init_bilan_eau
 
@@ -87 +91 @@
 endwhere
 
+archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel
+gr_line(:,:)=0
+do j=1,ny
+  do i=1,nx
+    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel)) then ! grounded with ice
+      if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) gr_line(i,j)=1
+      if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0..and..not.flot_marais(i+1,j)) gr_line(i,j)=1
+      if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0..and..not.flot_marais(i,j-1)) gr_line(i,j)=1
+      if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0..and..not.flot_marais(i,j+1)) gr_line(i,j)=1
+    endif
+  enddo
+enddo
+
+where (gr_line(:,:).eq.1)
+!~   grline_dtt(:,:)= (((uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))**2+ &
+!~                     (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))**2)**0.5)*0.5 &
+!~                     *H(:,:) + grline_dtt(:,:)
+  grline_dtt(:,:)= (abs(uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))/2.*dy + &
+                    abs(uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))/2.*dx)  &
+                    *H(:,:) + grline_dtt(:,:)             
+endwhere
+
+
 if (isynchro.eq.1) then
 ! on raisonne en bilan annuel pour simplifier :

trunk/SOURCES/initial-0.3.f90

@@ -30 +30 @@
   use sorties_ncdf_grisli
   use util_recovery
+  use bilan_flux_mod
 
   !-------------------------------------------------------------------------------------
@@ -135 +136 @@
   call init_calving
 
+  call init_bilan_flux
 !
 !-------------------------------------------------------------------------------------

trunk/SOURCES/steps_time_loop.f90

@@ -22 +22 @@
   use icetempmod
   use diagno_mod
-  use bilan_eau_mod 
+  use bilan_eau_mod
+  use bilan_flux_mod
   !  use track_debug 
 
@@ -64 +65 @@
 
      call bilan_eau
+     call bilan_flux_output
      if (isynchro.eq.1) then
          call shortoutput
@@ -72 +74 @@
          ablbord_dtt(:,:)=0.
          diff_H_2D(:,:)=0.
+         grline_dtt(:,:)=0.
       endif