Changeset 146 for branches/iLoveclim/SOURCES/Ant40_files/output_anta40_mod-0.4.f90

Timestamp:

10/20/17 09:31:39 (7 years ago)

Author:

aquiquet

Message:

Grisli-iLoveclim branch: merged to trunk at revision 145

Location:

branches/iLoveclim

Files:

• . (modified) (1 prop)
• SOURCES/Ant40_files/output_anta40_mod-0.4.f90 (modified) (8 diffs)

branches/iLoveclim/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)
+    read(568,'(a100)') standard_name(k)
+    read(568,'(a100)') long_name(k)
+    read(568,'(a100)') units(k)
+    read(568,*)
+  enddo
+  close(568)
+! 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)
+
+! definition des dimension :
+  status=nf90_def_dim(ncid, name="time", len=NF90_UNLIMITED, dimid=timeDimID)
+  if (status /= nf90_noerr) call handle_err(status)
+  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)
+  status=nf90_put_att(ncid, timeVarID,"units", "years since 2007-01-01 00:00:00")
+  if (status /= nf90_noerr) call handle_err(status)
+
+! 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))
+    if (status /= nf90_noerr) call handle_err(status)
+    status=nf90_put_att(ncid, varID(k), "standard_name", 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)))
+    if (status /= nf90_noerr) call handle_err(status)
+    status=nf90_put_att(ncid, varID(k), "units", trim(units(k)))
+    if (status /= nf90_noerr) call handle_err(status)
+  enddo
+
+! 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 +147 @@
 !------------------
 real ::  smax
+
       vol=0. 
       np=0
+      nflot=0
       hmax=0. 
       smax=0.
@@ -74 +169 @@
       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 +224 @@
         hmean=vol/np 
         vol=vol*dx*dy 
-        volf=volf*dx*dy
+        volf=volf*dx*dy*ice_density
         bmean=bmean/np 
         accmean=accmean/np 
@@ -130 +238 @@
         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/secyear
+      tendlibmassbf=tendlibmassbf*dx*dy*ice_density/secyear
+      tendlicalvf=tendlicalvf*dx*dy*ice_density/secyear
+      tendligroundf=tendligroundf*ice_density/secyear
+      
+      
       bdotmean=bdotmean/nx/ny  
 
@@ -154 +270 @@
 !940   format('%%%% ',a,'   time=',f8.0,' %%%%')
 
+
+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
+  
+  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