Index: /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 (revision 10957)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 (revision 10958)
@@ -35,4 +35,5 @@
USE phycst ! physical constants
USE in_out_manager ! I/O manager
+ USE iom ! I/0 library
USE daymod ! calendar
USE dianam ! build name of file
@@ -65,15 +66,22 @@
!! * Shared module variables
LOGICAL, PUBLIC, PARAMETER :: lk_diadct = .TRUE. !: model-data diagnostics flag
+ LOGICAL, PUBLIC :: ln_dct_calc_noos_25h !: Calcuate noos 25 h means
+ LOGICAL, PUBLIC :: ln_dct_calc_noos_hr !: Calcuate noos hourly means
+ ! JT
+ LOGICAL, PUBLIC :: ln_dct_iom_cont !: Use IOM Output?
+ LOGICAL, PUBLIC :: ln_dct_ascii !: Output ascii or binary
+ LOGICAL, PUBLIC :: ln_dct_h !: Output hourly instantaneous or mean values
+ ! JT
!! * Module variables
- INTEGER :: nn_dct = 1 ! Frequency of computation
- INTEGER :: nn_dctwri = 1 ! Frequency of output
- INTEGER :: nn_secdebug = 0 ! Number of the section to debug
- INTEGER :: nn_dct_h = 1 ! Frequency of computation for NOOS hourly files
- INTEGER :: nn_dctwri_h = 1 ! Frequency of output for NOOS hourly files
+ INTEGER :: nn_dct ! Frequency of computation
+ INTEGER :: nn_dctwri ! Frequency of output
+ INTEGER :: nn_secdebug ! Number of the section to debug
+ INTEGER :: nn_dct_h ! Frequency of computation for NOOS hourly files
+ INTEGER :: nn_dctwri_h ! Frequency of output for NOOS hourly files
INTEGER, PARAMETER :: nb_class_max = 12 ! maximum number of classes, i.e. depth levels or density classes
- INTEGER, PARAMETER :: nb_sec_max = 30 ! maximum number of sections
- INTEGER, PARAMETER :: nb_point_max = 300 ! maximum number of points in a single section
+ INTEGER, PARAMETER :: nb_sec_max = 100 ! maximum number of sections
+ INTEGER, PARAMETER :: nb_point_max = 375 ! maximum number of points in a single section
INTEGER, PARAMETER :: nb_type_class = 14 ! types of calculations, i.e. pos transport, neg transport, heat transport, salt transport
INTEGER, PARAMETER :: nb_3d_vars = 5
@@ -134,5 +142,5 @@
ALLOCATE(transports_3d_h(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(3) )
ALLOCATE(transports_2d_h(nb_2d_vars,nb_sec_max,nb_point_max) , STAT=ierr(4) )
- ALLOCATE(z_hr_output(nb_sec_max,168,nb_class_max) , STAT=ierr(5) ) ! 168 = 24 * 7days
+ ALLOCATE(z_hr_output(nb_sec_max,3,nb_class_max) , STAT=ierr(5) )
diadct_alloc = MAXVAL( ierr )
@@ -149,6 +157,7 @@
!!
!!---------------------------------------------------------------------
- NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug
- INTEGER :: ios ! Local integer output status for namelist read
+ NAMELIST/namdct/nn_dct,ln_dct_h,nn_dctwri,ln_dct_ascii,nn_secdebug,ln_dct_calc_noos_25h,ln_dct_calc_noos_hr,ln_dct_iom_cont
+ INTEGER :: ios,jsec ! Local integer output status for namelist read
+ CHARACTER(len=3) :: jsec_str ! name of the jsec
IF( nn_timing == 1 ) CALL timing_start('dia_dct_init')
@@ -164,7 +173,31 @@
IF( ln_NOOS ) THEN
- nn_dct=3600./rdt ! hard coded for NOOS transects, to give 25 hour means
+
+ !Do calculation for daily, 25hourly mean every hour
+ nn_dct=3600./rdt ! hard coded for NOOS transects, to give 25 hour means from hourly instantaneous values
+
+ !write out daily, 25hourly mean every day
nn_dctwri=86400./rdt
- nn_dct_h=1 ! hard coded for NOOS transects, to give hourly data
+
+
+ !nn_dct_h=1 ! hard coded for NOOS transects, to give hourly data
+ ! If you want hourly instantaneous values, you only do the calculation every 12 timesteps (if rdt = 300)
+ ! and output it every 12 time steps. For this, you set the ln_dct_h to be True, and it calcuates it automatically
+ ! if you want hourly mean values, set ln_dct_h to be False, and it will do the calculate every time step.
+ !
+ !SELECT CASE( ln_dct_h )
+ ! CASE(.TRUE.)
+ ! nn_dct_h=3600./rdt
+ ! CASE(.FALSE.)
+ ! nn_dct_h=1
+ !END SELECT
+
+ IF ( ln_dct_h ) THEN
+ nn_dct_h=3600./rdt
+ ELSE
+ nn_dct_h=1.
+ ENDIF
+
+ !JT write out hourly calculation every hour
nn_dctwri_h=3600./rdt
ENDIF
@@ -175,7 +208,14 @@
WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
IF( ln_NOOS ) THEN
+ WRITE(numout,*) " Calculate NOOS hourly output: ln_dct_calc_noos_hr = ",ln_dct_calc_noos_hr
+ WRITE(numout,*) " Calculate NOOS 25 hour mean output: ln_dct_calc_noos_hr = ",ln_dct_calc_noos_25h
+ WRITE(numout,*) " Use IOM Output: ln_dct_iom_cont = ",ln_dct_iom_cont
+ WRITE(numout,*) " Output in ASCII (True) or Binary (False): ln_dct_ascii = ",ln_dct_ascii
+ WRITE(numout,*) " Frequency of hourly computation - instantaneous (TRUE) or hourly mean (FALSE): ln_dct_h = ",ln_dct_h
+
WRITE(numout,*) " Frequency of computation hard coded to be every hour: nn_dct = ",nn_dct
WRITE(numout,*) " Frequency of write hard coded to average 25 instantaneous hour values: nn_dctwri = ",nn_dctwri
- WRITE(numout,*) " Frequency of hourly computation hard coded to be every timestep: nn_dct_h = ",nn_dct_h
+ WRITE(numout,*) " Frequency of hourly computation (timestep) : nn_dct_h = ",nn_dct_h
+ WRITE(numout,*) " Frequency of hourly computation Not hard coded to be every timestep, or : nn_dct_h = ",nn_dct_h
WRITE(numout,*) " Frequency of hourly write hard coded to every hour: nn_dctwri_h = ",nn_dctwri_h
ELSE
@@ -195,13 +235,22 @@
ENDIF
-
- !Read section_ijglobal.diadct
- CALL readsec
+
+
+ IF ( ln_NOOS ) THEN
+ IF ( ln_dct_calc_noos_25h .or. ln_dct_calc_noos_hr ) CALL readsec
+ ENDIF
!open output file
- IF( lwm ) THEN
+ IF( lwp ) THEN
IF( ln_NOOS ) THEN
- CALL ctl_opn( numdct_NOOS ,'NOOS_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
- CALL ctl_opn( numdct_NOOS_h,'NOOS_transport_h', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
+ WRITE(numout,*) "diadct_init: Open output files. ASCII? ",ln_dct_ascii
+ WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
+ IF ( ln_dct_ascii ) THEN
+ if ( ln_dct_calc_noos_25h ) CALL ctl_opn( numdct_NOOS ,'NOOS_transport' , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. )
+ if ( ln_dct_calc_noos_hr ) CALL ctl_opn( numdct_NOOS_h,'NOOS_transport_h', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. )
+ ELSE
+ if ( ln_dct_calc_noos_25h ) CALL ctl_opn( numdct_NOOS ,'NOOS_transport_bin' , 'REPLACE', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. )
+ if ( ln_dct_calc_noos_hr ) CALL ctl_opn( numdct_NOOS_h,'NOOS_transport_bin_h', 'REPLACE', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. )
+ ENDIF
ELSE
CALL ctl_opn( numdct_vol, 'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
@@ -212,11 +261,67 @@
! Initialise arrays to zero
- transports_3d(:,:,:,:) =0.0
- transports_2d(:,:,:) =0.0
- transports_3d_h(:,:,:,:)=0._wp
- transports_2d_h(:,:,:) =0._wp
- z_hr_output(:,:,:) =0._wp
+ transports_3d(:,:,:,:) =0._wp
+ transports_2d(:,:,:) =0._wp
+ transports_3d_h(:,:,:,:) =0._wp
+ transports_2d_h(:,:,:) =0._wp
+ z_hr_output(:,:,:) =0._wp
IF( nn_timing == 1 ) CALL timing_stop('dia_dct_init')
+
+ IF (ln_dct_iom_cont) THEN
+ IF( lwp ) THEN
+ WRITE(numout,*) " "
+ WRITE(numout,*) "diadct_init: using xios iom_put for output: field_def.xml and iodef.xml code"
+ WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
+ WRITE(numout,*) ""
+ WRITE(numout,*) " field_def.xml"
+ WRITE(numout,*) " ~~~~~~~~~~~~~"
+ WRITE(numout,*) ""
+ WRITE(numout,*) ""
+
+ WRITE(numout,*) ' '
+
+ DO jsec=1,nb_sec
+ WRITE (jsec_str, "(I3.3)") jsec
+
+ WRITE(numout,*) ' '
+ WRITE(numout,*) ' '
+ WRITE(numout,*) ' '
+
+ ENDDO
+
+ WRITE(numout,*) ' '
+
+ WRITE(numout,*) ""
+ WRITE(numout,*) ""
+ WRITE(numout,*) " iodef.xml"
+ WRITE(numout,*) " ~~~~~~~~~"
+ WRITE(numout,*) ""
+ WRITE(numout,*) ""
+
+ WRITE(numout,*) ' '
+ WRITE(numout,*) ""
+ WRITE(numout,*) ' '
+ DO jsec=1,nb_sec
+ WRITE (jsec_str, "(I3.3)") jsec
+
+ WRITE(numout,*) ' '
+ WRITE(numout,*) ' '
+ WRITE(numout,*) ' '
+
+ ENDDO
+ WRITE(numout,*) ' '
+ WRITE(numout,*) ""
+ WRITE(numout,*) ' '
+
+ WRITE(numout,*) ""
+ WRITE(numout,*) ""
+ WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
+ WRITE(numout,*) ""
+
+ ENDIF
+ ENDIF
+
+
!
END SUBROUTINE dia_dct_init
@@ -231,5 +336,5 @@
!! Method :: All arrays initialised to zero in dct_init
!! Each nn_dct time step call subroutine 'transports' for
- !! each section to sum the transports over each grid cell.
+ !! each section to sum the transports.
!! Each nn_dctwri time step:
!! Divide the arrays by the number of summations to gain
@@ -271,116 +376,148 @@
WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~~~~~"
WRITE(numout,*) "nb_sec = ",nb_sec
+ WRITE(numout,*) "nn_dct = ",nn_dct
+ WRITE(numout,*) "ln_NOOS = ",ln_NOOS
+ WRITE(numout,*) "nb_sec = ",nb_sec
+ WRITE(numout,*) "nb_sec_max = ",nb_sec_max
+ WRITE(numout,*) "nb_type_class = ",nb_type_class
+ WRITE(numout,*) "nb_class_max = ",nb_class_max
ENDIF
-
- ! Compute transport and write only at nn_dctwri
- IF( MOD(kt,nn_dct)==0 .or. & ! compute transport every nn_dct time steps
- (ln_NOOS .and. kt==nn_it000 ) ) THEN ! also include first time step when calculating NOOS 25 hour averages
-
- DO jsec=1,nb_sec
-
- !debug this section computing ?
- IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE.
-
- !Compute transport through section
- CALL transport(secs(jsec),lldebug,jsec)
-
- ENDDO
-
- IF( MOD(kt,nn_dctwri)==0 )THEN
-
- IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)" diadct: average transports and write at kt = ",kt
-
- !! divide arrays by nn_dctwri/nn_dct to obtain average
- transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct)
- transports_2d(:,:,:) =transports_2d(:,:,:) /(nn_dctwri/nn_dct)
-
- ! Sum over each class
- DO jsec=1,nb_sec
- CALL dia_dct_sum(secs(jsec),jsec)
- ENDDO
-
- !Sum on all procs
- IF( lk_mpp )THEN
- zsum(:,:,:)=0.0_wp
- ish(1) = nb_sec_max*nb_type_class*nb_class_max
- ish2 = (/nb_sec_max,nb_type_class,nb_class_max/)
- DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO
- zwork(:)= RESHAPE(zsum(:,:,:), ish )
- CALL mpp_sum(zwork, ish(1))
- zsum(:,:,:)= RESHAPE(zwork,ish2)
- DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO
- ENDIF
-
- !Write the transport
- DO jsec=1,nb_sec
-
- IF( lwm .and. .not. ln_NOOS )CALL dia_dct_wri(kt,jsec,secs(jsec))
- IF( lwm .and. ln_NOOS )CALL dia_dct_wri_NOOS(kt,jsec,secs(jsec)) ! use NOOS specific formatting
+
+ IF ( ln_dct_calc_noos_25h ) THEN
+
+ ! Compute transport and write only at nn_dctwri
+ IF ( MOD(kt,nn_dct)==0 .or. & ! compute transport every nn_dct time steps
+ (ln_NOOS .and. kt==nn_it000 ) ) THEN ! also include first time step when calculating NOOS 25 hour averages
- !nullify transports values after writing
- transports_3d(:,jsec,:,:)=0.
- transports_2d(:,jsec,: )=0.
- secs(jsec)%transport(:,:)=0.
- IF ( ln_NOOS ) CALL transport(secs(jsec),lldebug,jsec) ! reinitialise for next 25 hour instantaneous average (overlapping values)
-
- ENDDO
-
- ENDIF
-
+
+
+ DO jsec=1,nb_sec
+
+ lldebug=.FALSE.
+ IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE.
+
+ !Compute transport through section
+ CALL transport(secs(jsec),lldebug,jsec)
+
+
+ ENDDO
+
+ IF( MOD(kt,nn_dctwri)==0 )THEN
+
+
+
+ IF( lwp .AND. kt==nit000+nn_dctwri-1 ) WRITE(numout,*)" diadct: average and write at kt = ",kt
+
+
+ ! Not 24 values, but 25! divide by ((nn_dctwri/nn_dct) +1)
+ !! divide arrays by nn_dctwri/nn_dct to obtain average
+ transports_3d(:,:,:,:)= transports_3d(:,:,:,:)/((nn_dctwri/nn_dct)+1.)
+ transports_2d(:,:,:) = transports_2d(:,:,:) /((nn_dctwri/nn_dct)+1.)
+
+ ! Sum over each class
+ DO jsec=1,nb_sec
+ CALL dia_dct_sum(secs(jsec),jsec)
+ ENDDO
+
+ !Sum on all procs
+ IF( lk_mpp )THEN
+ zsum(:,:,:)=0.0_wp
+ ish(1) = nb_sec_max*nb_type_class*nb_class_max
+ ish2 = (/nb_sec_max,nb_type_class,nb_class_max/)
+ DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO
+ zwork(:)= RESHAPE(zsum(:,:,:), ish )
+ CALL mpp_sum(zwork, ish(1))
+ zsum(:,:,:)= RESHAPE(zwork,ish2)
+ DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO
+ ENDIF
+
+ !Write the transport
+ DO jsec=1,nb_sec
+
+ IF( lwp .and. .not. ln_NOOS )CALL dia_dct_wri(kt,jsec,secs(jsec))
+ !IF( lwp .and. ln_NOOS )CALL dia_dct_wri_NOOS(kt,jsec,secs(jsec)) ! use NOOS specific formatting
+ IF( ln_NOOS )CALL dia_dct_wri_NOOS(kt,jsec,secs(jsec)) ! use NOOS specific formatting
+
+
+ !nullify transports values after writing
+ transports_3d(:,jsec,:,:)=0.0
+ transports_2d(:,jsec,: )=0.0
+ secs(jsec)%transport(:,:)=0.
+
+
+ IF ( ln_NOOS ) CALL transport(secs(jsec),lldebug,jsec) ! reinitialise for next 25 hour instantaneous average (overlapping values)
+
+
+
+ ENDDO
+
+ ENDIF
+
+ ENDIF
+
+ ENDIF
+ IF ( ln_dct_calc_noos_hr ) THEN
+ IF ( MOD(kt,nn_dct_h)==0 ) THEN ! compute transport every nn_dct_h time steps
+
+ DO jsec=1,nb_sec
+
+ lldebug=.FALSE.
+ IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. kt==nit000+nn_dct_h-1 .AND. lwp ) lldebug=.TRUE.
+
+ !Compute transport through section
+ CALL transport_h(secs(jsec),lldebug,jsec)
+
+ ENDDO
+
+ IF( MOD(kt,nn_dctwri_h)==0 )THEN
+
+ IF( lwp .AND. kt==nit000+nn_dctwri_h-1 )WRITE(numout,*)" diadct: average and write hourly files at kt = ",kt
+
+ !! divide arrays by nn_dctwri/nn_dct to obtain average
+ !
+ ! JT - I think this is wrong. I think it is trying to sum over 25 hours, but only dividing by 24.
+ ! I think it might work for daily cycles, but not for monthly cycles,
+ !
+ transports_3d_h(:,:,:,:)=transports_3d_h(:,:,:,:)/(nn_dctwri_h/nn_dct_h)
+ transports_2d_h(:,:,:) =transports_2d_h(:,:,:) /(nn_dctwri_h/nn_dct_h)
+
+ ! Sum over each class
+ DO jsec=1,nb_sec
+ CALL dia_dct_sum_h(secs(jsec),jsec)
+ ENDDO
+
+ !Sum on all procs
+ IF( lk_mpp )THEN
+ ish(1) = nb_sec_max*nb_type_class*nb_class_max
+ ish2 = (/nb_sec_max,nb_type_class,nb_class_max/)
+ DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport_h(:,:) ; ENDDO
+ zwork(:)= RESHAPE(zsum(:,:,:), ish )
+ CALL mpp_sum(zwork, ish(1))
+ zsum(:,:,:)= RESHAPE(zwork,ish2)
+ DO jsec=1,nb_sec ; secs(jsec)%transport_h(:,:) = zsum(jsec,:,:) ; ENDDO
+ ENDIF
+
+ !Write the transport
+ DO jsec=1,nb_sec
+
+ IF( lwp .and. ln_NOOS ) THEN
+ CALL dia_dct_wri_NOOS_h(kt/nn_dctwri_h,jsec,secs(jsec)) ! use NOOS specific formatting
+ endif
+ !nullify transports values after writing
+ transports_3d_h(:,jsec,:,:)=0.0
+ transports_2d_h(:,jsec,:)=0.0
+ secs(jsec)%transport_h(:,:)=0.0
+
+ ! for hourly mean or hourly instantaneous, you don't initialise! start with zero!
+ !IF ( ln_NOOS ) CALL transport_h(secs(jsec),lldebug,jsec) ! reinitialise for next 25 hour instantaneous average (overlapping values)
+
+ ENDDO
+
+ ENDIF
+
+ ENDIF
+
ENDIF
-
- IF ( MOD(kt,nn_dct_h)==0 ) THEN ! compute transport every nn_dct_h time steps
-
- DO jsec=1,nb_sec
-
- !lldebug=.FALSE.
- IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. kt==nit000+nn_dct_h-1 .AND. lwp ) lldebug=.TRUE.
-
- !Compute transport through section
- CALL transport_h(secs(jsec),lldebug,jsec)
-
- ENDDO
-
- IF( MOD(kt,nn_dctwri_h)==0 )THEN
-
- IF( lwp .AND. kt==nit000+nn_dctwri_h-1 )WRITE(numout,*)" diadct: average and write hourly files at kt = ",kt
-
- !! divide arrays by nn_dctwri/nn_dct to obtain average
- transports_3d_h(:,:,:,:)=transports_3d_h(:,:,:,:)/(nn_dctwri_h/nn_dct_h)
- transports_2d_h(:,:,:) =transports_2d_h(:,:,:) /(nn_dctwri_h/nn_dct_h)
-
- ! Sum over each class
- DO jsec=1,nb_sec
- CALL dia_dct_sum_h(secs(jsec),jsec)
- ENDDO
-
- !Sum on all procs
- IF( lk_mpp )THEN
- ish(1) = nb_sec_max*nb_type_class*nb_class_max
- ish2 = (/nb_sec_max,nb_type_class,nb_class_max/)
- DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport_h(:,:) ; ENDDO
- zwork(:)= RESHAPE(zsum(:,:,:), ish )
- CALL mpp_sum(zwork, ish(1))
- zsum(:,:,:)= RESHAPE(zwork,ish2)
- DO jsec=1,nb_sec ; secs(jsec)%transport_h(:,:) = zsum(jsec,:,:) ; ENDDO
- ENDIF
-
- !Write the transport
- DO jsec=1,nb_sec
-
- IF( lwp .and. ln_NOOS )CALL dia_dct_wri_NOOS_h(kt/nn_dctwri_h,jsec,secs(jsec)) ! use NOOS specific formatting
-
- !nullify transports values after writing
- transports_3d_h(:,jsec,:,:)=0.0
- transports_2d_h(:,jsec,:)=0.0
- secs(jsec)%transport_h(:,:)=0.
- IF ( ln_NOOS ) CALL transport_h(secs(jsec),lldebug,jsec) ! reinitialise for next 25 hour instantaneous average (overlapping values)
-
- ENDDO
-
- ENDIF
-
- ENDIF
IF( lk_mpp )THEN
@@ -424,5 +561,20 @@
!open input file
!---------------
- CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
+ !write(numout,*) 'dct low-level pre open: little endian '
+ !OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='LITTLE_ENDIAN')
+
+ write(numout,*) 'dct low-level pre open: big endian :',nproc,narea
+ OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='BIG_ENDIAN')
+
+ !write(numout,*) 'dct low-level pre open: SWAP '
+ !OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='SWAP')
+
+ !write(numout,*) 'dct low-level pre open: NATIVE '
+ !OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='NATIVE')
+
+ READ(107) isec
+ CLOSE(107)
+
+ CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. )
!---------------
@@ -433,25 +585,36 @@
IF ( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) ) &
- & WRITE(numout,*)'debuging for section number: ',jsec
+ & WRITE(numout,*)'debugging for section number: ',jsec
!initialization
!---------------
secs(jsec)%name=''
- secs(jsec)%llstrpond = .FALSE. ; secs(jsec)%ll_ice_section = .FALSE.
- secs(jsec)%ll_date_line = .FALSE. ; secs(jsec)%nb_class = 0
- secs(jsec)%zsigi = 99._wp ; secs(jsec)%zsigp = 99._wp
- secs(jsec)%zsal = 99._wp ; secs(jsec)%ztem = 99._wp
- secs(jsec)%zlay = 99._wp
- secs(jsec)%transport = 0._wp ; secs(jsec)%nb_point = 0
- secs(jsec)%transport_h = 0._wp ; secs(jsec)%nb_point = 0
+ secs(jsec)%llstrpond = .FALSE.
+ secs(jsec)%ll_ice_section = .FALSE.
+ secs(jsec)%ll_date_line = .FALSE.
+ secs(jsec)%nb_class = 0
+ secs(jsec)%zsigi = 99._wp
+ secs(jsec)%zsigp = 99._wp
+ secs(jsec)%zsal = 99._wp
+ secs(jsec)%ztem = 99._wp
+ secs(jsec)%zlay = 99._wp
+ secs(jsec)%transport = 0._wp
+ secs(jsec)%transport_h = 0._wp
+ secs(jsec)%nb_point = 0
!read section's number / name / computing choices / classes / slopeSection / points number
!-----------------------------------------------------------------------------------------
- READ(numdct_in,iostat=iost)isec
- IF (iost .NE. 0 )EXIT !end of file
+
+ READ(numdct_in,iostat=iost) isec
+ IF (iost .NE. 0 ) then
+ write(numout,*) 'unable to read section_ijglobal.diadct. iost = ',iost
+ EXIT !end of file
+ ENDIF
+
+
WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec
- IF( jsec .NE. isec )CALL ctl_stop( cltmp )
-
- IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )WRITE(numout,*)"isec ",isec
+
+
+ IF( jsec .NE. isec ) CALL ctl_stop( cltmp )
READ(numdct_in)secs(jsec)%name
@@ -483,5 +646,5 @@
WRITE(numout,*) " Section name : ",TRIM(secs(jsec)%name)
- WRITE(numout,*) " Compute heat and salt transport ? ",secs(jsec)%llstrpond
+ WRITE(numout,*) " Compute temperature and salinity transports ? ",secs(jsec)%llstrpond
WRITE(numout,*) " Compute ice transport ? ",secs(jsec)%ll_ice_section
WRITE(numout,*) " Section crosses date-line ? ",secs(jsec)%ll_date_line
@@ -558,10 +721,10 @@
ENDIF
- IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
+ IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
DO jpt = 1,iptloc
iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
ENDDO
- ENDIF
+ ENDIF
!remove redundant points between processors
@@ -602,4 +765,6 @@
nb_sec = jsec-1 !number of section read in the file
+
+ IF( lwp ) WRITE(numout,*)'diadct: read sections: Finished readsec.'
CALL wrk_dealloc( nb_point_max, directemp )
@@ -703,6 +868,6 @@
!! loop on the level jk !!
!!
- !! Output :: Arrays containing the volume,density,heat,salt transports for each i
- !! point in a section, summed over each nn_dct.
+ !! ** Output: Arrays containing the volume,density,salinity,temperature etc
+ !! transports for each point in a section, summed over each nn_dct.
!!
!!-------------------------------------------------------------------------------------------
@@ -713,15 +878,15 @@
!! * Local variables
- INTEGER :: jk, jseg, jclass,jl, &!loop on level/segment/classes/ice categories
- isgnu, isgnv !
- REAL(wp) :: zumid, zvmid, &!U/V velocity on a cell segment
- zumid_ice, zvmid_ice, &!U/V ice velocity
- zTnorm !transport of velocity through one cell's sides
- REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep !temperature/salinity/potential density/ssh/depth at u/v point
+ INTEGER :: jk,jseg,jclass, &!loop on level/segment/classes
+ isgnu , isgnv !
+ REAL(wp):: zumid , zvmid , &!U/V velocity on a cell segment
+ zumid_ice , zvmid_ice , &!U/V ice velocity
+ zTnorm !transport of velocity through one cell's sides
+ REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
TYPE(POINT_SECTION) :: k
!!--------------------------------------------------------
- !!NIALL IF( ld_debug )WRITE(numout,*)' Compute transport'
+ IF( ld_debug )WRITE(numout,*)' Compute transport (jsec,sec%nb_point,sec%slopeSection) : ', jsec,sec%nb_point,sec%slopeSection
!---------------------------!
@@ -730,4 +895,18 @@
IF(sec%nb_point .NE. 0)THEN
+ !----------------------------------------------------------------------------------------------------
+ !----------------------------------------------------------------------------------------------------
+ !----------------------------------------------------------------------------------------------------
+ !
+ !
+ ! ! ! ! JT 1/09/2018 - changing convention. Always direction + is toward left hand of section
+ !
+ ! Making sign of the velocities used to calculate the volume transport a function of direction, not slopesection
+ ! (isgnu, isgnv)
+ !
+ ! They vary for each segment of the section.
+ !
+ !----------------------------------------------------------------------------------------------------
+ !----------------------------------------------------------------------------------------------------
!----------------------------------------------------------------------------------------------------
!Compute sign for velocities:
@@ -751,9 +930,4 @@
!
!----------------------------------------------------------------------------------------------------
- isgnu = 1
- IF( sec%slopeSection .GT. 0 ) THEN ; isgnv = -1
- ELSE ; isgnv = 1
- ENDIF
- IF( sec%slopeSection .GE. 9999. ) isgnv = 1
IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
@@ -763,31 +937,62 @@
!--------------------------------------!
DO jseg=1,MAX(sec%nb_point-1,0)
-
+
+
+ !Compute sign for velocities:
+
+ !isgnu = 1
+ !isgnv = 1
+ !
+ !changing sign of u and v is dependent on the direction of the section.
+ !isgnu = 1
+ !isgnv = 1
+ !SELECT CASE( sec%direction(jseg) )
+ !CASE(0) ; isgnv = -1
+ !CASE(3) ; isgnu = -1
+ !END SELECT
+
+
+ SELECT CASE( sec%direction(jseg) )
+ CASE(0)
+ isgnu = 1
+ isgnv = -1
+ CASE(1)
+ isgnu = 1
+ isgnv = 1
+ CASE(2)
+ isgnu = 1
+ isgnv = 1
+ CASE(3)
+ isgnu = -1
+ isgnv = 1
+ END SELECT
+
!-------------------------------------------------------------------------------------------
! Select the appropriate coordinate for computing the velocity of the segment
+ ! Corrected by JT 01/09/2018 (#)
!
! CASE(0) Case (2)
! ------- --------
! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j)
- ! F(i,j)----------V(i+1,j)-------F(i+1,j) |
- ! |
- ! |
- ! |
- ! Case (3) U(i,j)
- ! -------- |
- ! |
+ ! F(i,j)---------#V(i,j)-------F(i+1,j) |
+ ! --------> |
+ ! | |
+ ! | |
+ ! Case (3) | U(i,j)
+ ! -------- | |
+ ! V |
! listPoint(jseg+1) F(i,j+1) |
! | |
! | |
! | listPoint(jseg+1) F(i,j-1)
- ! |
- ! |
- ! U(i,j+1)
- ! | Case(1)
- ! | ------
+ ! ^ |
+ ! | |
+ ! | U(i,j+1)
+ ! | | Case(1)
+ ! | | ------
! |
! | listPoint(jseg+1) listPoint(jseg)
- ! | F(i-1,j)-----------V(i,j) -------f(jseg)
- ! listPoint(jseg) F(i,j)
+ ! | F(i-1,j)----------#V(i-1,j) ------#f(i,j)
+ ! listPoint(jseg) F(i,j) <-------
!
!-------------------------------------------------------------------------------------------
@@ -800,57 +1005,57 @@
END SELECT
- !---------------------------|
- ! LOOP ON THE LEVEL |
- !---------------------------|
- !Sum of the transport on the vertical
- DO jk=1,mbathy(k%I,k%J)
-
- ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
- SELECT CASE( sec%direction(jseg) )
- CASE(0,1)
- ztn = interp(k%I,k%J,jk,'V',0 )
- zsn = interp(k%I,k%J,jk,'V',1 )
- zrhop = interp(k%I,k%J,jk,'V',2 )
- zrhoi = interp(k%I,k%J,jk,'V',3 )
- zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1)
- CASE(2,3)
- ztn = interp(k%I,k%J,jk,'U',0 )
- zsn = interp(k%I,k%J,jk,'U',1 )
- zrhop = interp(k%I,k%J,jk,'U',2 )
- zrhoi = interp(k%I,k%J,jk,'U',3 )
- zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1)
- END SELECT
-
- zfsdep= fsdept(k%I,k%J,jk)
-
- !compute velocity with the correct direction
- SELECT CASE( sec%direction(jseg) )
- CASE(0,1)
- zumid=0.
- zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)
- CASE(2,3)
- zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)
- zvmid=0.
- END SELECT
-
- !zTnorm=transport through one cell;
- !velocity* cell's length * cell's thickness
- zTnorm=zumid*e2u(k%I,k%J)* fse3u(k%I,k%J,jk)+ &
- zvmid*e1v(k%I,k%J)* fse3v(k%I,k%J,jk)
+ !---------------------------|
+ ! LOOP ON THE LEVEL |
+ !---------------------------|
+ !Sum of the transport on the vertical
+ DO jk=1,mbathy(k%I,k%J)
+
+ ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
+ SELECT CASE( sec%direction(jseg) )
+ CASE(0,1)
+ ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,0) )
+ zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,1) )
+ zrhop = interp(k%I,k%J,jk,'V',2)
+ zrhoi = interp(k%I,k%J,jk,'V',3)
+ zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1)
+ CASE(2,3)
+ ztn = interp(k%I,k%J,jk,'U',0)
+ zsn = interp(k%I,k%J,jk,'U',1)
+ zrhop = interp(k%I,k%J,jk,'U',2)
+ zrhoi = interp(k%I,k%J,jk,'U',3)
+ zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1)
+ END SELECT
+
+ zfsdep= fsdept(k%I,k%J,jk)
+
+ !compute velocity with the correct direction
+ SELECT CASE( sec%direction(jseg) )
+ CASE(0,1)
+ zumid=0.
+ zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)
+ CASE(2,3)
+ zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)
+ zvmid=0.
+ END SELECT
+
+ !zTnorm=transport through one cell;
+ !velocity* cell's length * cell's thickness
+ zTnorm=zumid*e2u(k%I,k%J)* fse3u(k%I,k%J,jk)+ &
+ zvmid*e1v(k%I,k%J)* fse3v(k%I,k%J,jk)
#if ! defined key_vvl
- !add transport due to free surface
- IF( jk==1 )THEN
- zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &
- zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
- ENDIF
+ !add transport due to free surface
+ IF( jk==1 )THEN
+ zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &
+ zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
+ ENDIF
#endif
!COMPUTE TRANSPORT
-
- transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm
-
- IF ( sec%llstrpond ) THEN
- transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk) + zTnorm * ztn * zrhop * rcp
- transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk) + zTnorm * zsn * zrhop * 0.001
+
+ transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm
+
+ IF ( sec%llstrpond ) THEN
+ transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk) + zTnorm * zrhoi
+ transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk) + zTnorm * zrhop
transports_3d(4,jsec,jseg,jk) = transports_3d(4,jsec,jseg,jk) + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp
transports_3d(5,jsec,jseg,jk) = transports_3d(5,jsec,jseg,jk) + zTnorm * 0.001 * zsn * 1026._wp
@@ -880,24 +1085,13 @@
zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
-
-#if defined key_lim2
- transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)* &
- (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
- *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + &
- hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))
- transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)* &
- (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))
-#endif
-#if defined key_lim3
- DO jl=1,jpl
- transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)* &
- a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) * &
- ( ht_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) + &
- ht_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) )
-
- transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)* &
- a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl)
- ENDDO
-#endif
+
+ transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)* &
+ (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
+ *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + &
+ hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
+ +zice_vol_pos
+ transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)* &
+ (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
+ +zice_surf_pos
ENDIF !end of ice case
@@ -957,4 +1151,18 @@
!----------------------------------------------------------------------------------------------------
+ !----------------------------------------------------------------------------------------------------
+ !----------------------------------------------------------------------------------------------------
+ !
+ !
+ ! ! ! ! JT 1/09/2018 - changing convention. Always direction + is toward left hand of section
+ !
+ ! Making sign of the velocities used to calculate the volume transport a function of direction, not slopesection
+ ! (isgnu, isgnv)
+ !
+ ! They vary for each segment of the section.
+ !
+ !----------------------------------------------------------------------------------------------------
+ !----------------------------------------------------------------------------------------------------
+ !----------------------------------------------------------------------------------------------------
!Compute sign for velocities:
!
@@ -977,8 +1185,4 @@
!
!----------------------------------------------------------------------------------------------------
- isgnu = 1
- IF( sec%slopeSection .GT. 0 ) THEN ; isgnv = -1
- ELSE ; isgnv = 1
- ENDIF
IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv
@@ -988,31 +1192,62 @@
!--------------------------------------!
DO jseg=1,MAX(sec%nb_point-1,0)
-
+
+
+ !Compute sign for velocities:
+
+ !isgnu = 1
+ !isgnv = 1
+ !
+ ! changing sign of u and v is dependent on the direction of the section.
+ !isgnu = 1
+ !isgnv = 1
+ !SELECT CASE( sec%direction(jseg) )
+ !CASE(0) ; isgnv = -1
+ !CASE(3) ; isgnu = -1
+ !END SELECT
+
+
+ SELECT CASE( sec%direction(jseg) )
+ CASE(0)
+ isgnu = 1
+ isgnv = -1
+ CASE(1)
+ isgnu = 1
+ isgnv = 1
+ CASE(2)
+ isgnu = 1
+ isgnv = 1
+ CASE(3)
+ isgnu = -1
+ isgnv = 1
+ END SELECT
+
!-------------------------------------------------------------------------------------------
! Select the appropriate coordinate for computing the velocity of the segment
+ ! Corrected by JT 01/09/2018 (#)
!
! CASE(0) Case (2)
! ------- --------
! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j)
- ! F(i,j)----------V(i+1,j)-------F(i+1,j) |
- ! |
- ! |
- ! |
- ! Case (3) U(i,j)
- ! -------- |
- ! |
+ ! F(i,j)---------#V(i,j)-------F(i+1,j) |
+ ! --------> |
+ ! | |
+ ! | |
+ ! Case (3) | U(i,j)
+ ! -------- | |
+ ! V |
! listPoint(jseg+1) F(i,j+1) |
! | |
! | |
! | listPoint(jseg+1) F(i,j-1)
- ! |
- ! |
- ! U(i,j+1)
- ! | Case(1)
- ! | ------
+ ! ^ |
+ ! | |
+ ! | U(i,j+1)
+ ! | | Case(1)
+ ! | | ------
! |
! | listPoint(jseg+1) listPoint(jseg)
- ! | F(i-1,j)-----------V(i,j) -------f(jseg)
- ! listPoint(jseg) F(i,j)
+ ! | F(i-1,j)----------#V(i-1,j) ------#f(i,j)
+ ! listPoint(jseg) F(i,j) <-------
!
!-------------------------------------------------------------------------------------------
@@ -1031,8 +1266,8 @@
DO jk=1,mbathy(k%I,k%J)
- ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point
+ ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
SELECT CASE( sec%direction(jseg) )
CASE(0,1)
- ztn = interp(k%I,k%J,jk,'V',0 )
+ ztn = interp(k%I,k%J,jk,'V',0)
zsn = interp(k%I,k%J,jk,'V',1)
zrhop = interp(k%I,k%J,jk,'V',2)
@@ -1151,215 +1386,5 @@
REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
!!-------------------------------------------------------------
-
- !! Sum the relevant segments to obtain values for each class
- IF(sec%nb_point .NE. 0)THEN
-
- !--------------------------------------!
- ! LOOP ON THE SEGMENT BETWEEN 2 NODES !
- !--------------------------------------!
- DO jseg=1,MAX(sec%nb_point-1,0)
-
- !-------------------------------------------------------------------------------------------
- ! Select the appropriate coordinate for computing the velocity of the segment
- !
- ! CASE(0) Case (2)
- ! ------- --------
- ! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j)
- ! F(i,j)----------V(i+1,j)-------F(i+1,j) |
- ! |
- ! |
- ! |
- ! Case (3) U(i,j)
- ! -------- |
- ! |
- ! listPoint(jseg+1) F(i,j+1) |
- ! | |
- ! | |
- ! | listPoint(jseg+1) F(i,j-1)
- ! |
- ! |
- ! U(i,j+1)
- ! | Case(1)
- ! | ------
- ! |
- ! | listPoint(jseg+1) listPoint(jseg)
- ! | F(i-1,j)-----------V(i,j) -------f(jseg)
- ! listPoint(jseg) F(i,j)
- !
- !-------------------------------------------------------------------------------------------
-
- SELECT CASE( sec%direction(jseg) )
- CASE(0) ; k = sec%listPoint(jseg)
- CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
- CASE(2) ; k = sec%listPoint(jseg)
- CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
- END SELECT
-
- !---------------------------|
- ! LOOP ON THE LEVEL |
- !---------------------------|
- !Sum of the transport on the vertical
- DO jk=1,mbathy(k%I,k%J)
-
- ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
- SELECT CASE( sec%direction(jseg) )
- CASE(0,1)
- ztn = interp(k%I,k%J,jk,'V',0 )
- zsn = interp(k%I,k%J,jk,'V',1 )
- zrhop = interp(k%I,k%J,jk,'V',2)
- zrhoi = interp(k%I,k%J,jk,'V',3)
-
- CASE(2,3)
- ztn = interp(k%I,k%J,jk,'U',0 )
- zsn = interp(k%I,k%J,jk,'U',1 )
- zrhop = interp(k%I,k%J,jk,'U',2 )
- zrhoi = interp(k%I,k%J,jk,'U',3 )
- zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1)
- END SELECT
-
- zfsdep= fsdept(k%I,k%J,jk)
-
- !-------------------------------
- ! LOOP ON THE DENSITY CLASSES |
- !-------------------------------
- !The computation is made for each density/temperature/salinity/depth class
- DO jclass=1,MAX(1,sec%nb_class-1)
-
- !----------------------------------------------!
- !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
- !----------------------------------------------!
-
- IF ( ( &
- ((( zrhop .GE. (sec%zsigp(jclass)+1000. )) .AND. &
- ( zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR. &
- ( sec%zsigp(jclass) .EQ. 99.)) .AND. &
-
- ((( zrhoi .GE. (sec%zsigi(jclass) + 1000. )) .AND. &
- ( zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR. &
- ( sec%zsigi(jclass) .EQ. 99.)) .AND. &
-
- ((( zsn .GT. sec%zsal(jclass)) .AND. &
- ( zsn .LE. sec%zsal(jclass+1))) .OR. &
- ( sec%zsal(jclass) .EQ. 99.)) .AND. &
-
- ((( ztn .GE. sec%ztem(jclass)) .AND. &
- ( ztn .LE. sec%ztem(jclass+1))) .OR. &
- ( sec%ztem(jclass) .EQ.99.)) .AND. &
-
- ((( zfsdep .GE. sec%zlay(jclass)) .AND. &
- ( zfsdep .LE. sec%zlay(jclass+1))) .OR. &
- ( sec%zlay(jclass) .EQ. 99. )) &
- )) THEN
-
- !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
- !----------------------------------------------------------------------------
- IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN
- sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6
- ELSE
- sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6
- ENDIF
- IF( sec%llstrpond )THEN
-
- IF ( transports_3d(2,jsec,jseg,jk) .GE. 0.0 ) THEN
- sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)
- ELSE
- sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)
- ENDIF
-
- IF ( transports_3d(3,jsec,jseg,jk) .GE. 0.0 ) THEN
- sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)
- ELSE
- sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)
- ENDIF
-
- IF ( transports_3d(4,jsec,jseg,jk) .GE. 0.0 ) THEN
- sec%transport(7,jclass) = sec%transport(7,jclass)+transports_3d(4,jsec,jseg,jk)
- ELSE
- sec%transport(8,jclass) = sec%transport(8,jclass)+transports_3d(4,jsec,jseg,jk)
- ENDIF
-
- IF ( transports_3d(5,jsec,jseg,jk) .GE. 0.0 ) THEN
- sec%transport( 9,jclass) = sec%transport( 9,jclass)+transports_3d(5,jsec,jseg,jk)
- ELSE
- sec%transport(10,jclass) = sec%transport(10,jclass)+transports_3d(5,jsec,jseg,jk)
- ENDIF
-
- ELSE
- sec%transport( 3,jclass) = 0._wp
- sec%transport( 4,jclass) = 0._wp
- sec%transport( 5,jclass) = 0._wp
- sec%transport( 6,jclass) = 0._wp
- sec%transport( 7,jclass) = 0._wp
- sec%transport( 8,jclass) = 0._wp
- sec%transport( 9,jclass) = 0._wp
- sec%transport(10,jclass) = 0._wp
- ENDIF
-
- ENDIF ! end of test if point is in class
-
- ENDDO ! end of loop on the classes
-
- ENDDO ! loop over jk
-
-#if defined key_lim2 || defined key_lim3
-
- !ICE CASE
- IF( sec%ll_ice_section )THEN
-
- IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN
- sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg)*1.E-6
- ELSE
- sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg)*1.E-6
- ENDIF
-
- IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN
- sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg)*1.E-6
- ELSE
- sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg)*1.E-6
- ENDIF
-
- ENDIF !end of ice case
-#endif
-
- ENDDO !end of loop on the segment
-
- ELSE !if sec%nb_point =0
- sec%transport(1:2,:)=0.
- IF (sec%llstrpond) sec%transport(3:10,:)=0.
- IF (sec%ll_ice_section) sec%transport(11:14,:)=0.
- ENDIF !end of sec%nb_point =0 case
-
- END SUBROUTINE dia_dct_sum
-
- SUBROUTINE dia_dct_sum_h(sec,jsec)
- !!-------------------------------------------------------------
- !! Exactly as dia_dct_sum but for hourly files containing data summed at each time step
- !!
- !! Purpose: Average the transport over nn_dctwri time steps
- !! and sum over the density/salinity/temperature/depth classes
- !!
- !! Method:
- !! Sum over relevant grid cells to obtain values
- !! for each
- !! There are several loops:
- !! loop on the segment between 2 nodes
- !! loop on the level jk
- !! loop on the density/temperature/salinity/level classes
- !! test on the density/temperature/salinity/level
- !!
- !! ** Method :Transport through a given section is equal to the sum of transports
- !! computed on each proc.
- !! On each proc,transport is equal to the sum of transport computed through
- !! segments linking each point of sec%listPoint with the next one.
- !!
- !!-------------------------------------------------------------
- !! * arguments
- TYPE(SECTION),INTENT(INOUT) :: sec
- INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section
-
- TYPE(POINT_SECTION) :: k
- INTEGER :: jk,jseg,jclass !loop on level/segment/classes
- REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
- !!-------------------------------------------------------------
+
!! Sum the relevant segments to obtain values for each class
@@ -1416,14 +1441,229 @@
SELECT CASE( sec%direction(jseg) )
CASE(0,1)
- ztn = interp(k%I,k%J,jk,'V',0 )
- zsn = interp(k%I,k%J,jk,'V',1 )
- zrhop = interp(k%I,k%J,jk,'V',2 )
- zrhoi = interp(k%I,k%J,jk,'V',3 )
+ ztn = interp(k%I,k%J,jk,'V',0)
+ zsn = interp(k%I,k%J,jk,'V',1)
+ zrhop = interp(k%I,k%J,jk,'V',2)
+ zrhoi = interp(k%I,k%J,jk,'V',3)
zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1)
CASE(2,3)
- ztn = interp(k%I,k%J,jk,'U',0 )
- zsn = interp(k%I,k%J,jk,'U',1 )
- zrhop = interp(k%I,k%J,jk,'U',2 )
- zrhoi = interp(k%I,k%J,jk,'U',3 )
+ ztn = interp(k%I,k%J,jk,'U',0)
+ zsn = interp(k%I,k%J,jk,'U',1)
+ zrhop = interp(k%I,k%J,jk,'U',2)
+ zrhoi = interp(k%I,k%J,jk,'U',3)
+ zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1)
+ END SELECT
+
+ zfsdep= fsdept(k%I,k%J,jk)
+
+ !-------------------------------
+ ! LOOP ON THE DENSITY CLASSES |
+ !-------------------------------
+ !The computation is made for each density/temperature/salinity/depth class
+ DO jclass=1,MAX(1,sec%nb_class-1)
+
+ !----------------------------------------------!
+ !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
+ !----------------------------------------------!
+
+ IF ( ( &
+ ((( zrhop .GE. (sec%zsigp(jclass)+1000. )) .AND. &
+ ( zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR. &
+ ( sec%zsigp(jclass) .EQ. 99.)) .AND. &
+
+ ((( zrhoi .GE. (sec%zsigi(jclass) + 1000. )) .AND. &
+ ( zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR. &
+ ( sec%zsigi(jclass) .EQ. 99.)) .AND. &
+
+ ((( zsn .GT. sec%zsal(jclass)) .AND. &
+ ( zsn .LE. sec%zsal(jclass+1))) .OR. &
+ ( sec%zsal(jclass) .EQ. 99.)) .AND. &
+
+ ((( ztn .GE. sec%ztem(jclass)) .AND. &
+ ( ztn .LE. sec%ztem(jclass+1))) .OR. &
+ ( sec%ztem(jclass) .EQ.99.)) .AND. &
+
+ ((( zfsdep .GE. sec%zlay(jclass)) .AND. &
+ ( zfsdep .LE. sec%zlay(jclass+1))) .OR. &
+ ( sec%zlay(jclass) .EQ. 99. )) &
+ )) THEN
+
+ !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
+ !----------------------------------------------------------------------------
+ IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN
+ sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)
+ ELSE
+ sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)
+ ENDIF
+ IF( sec%llstrpond )THEN
+
+ IF( transports_3d(1,jsec,jseg,jk) .NE. 0._wp ) THEN
+
+ IF (transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN
+ sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
+ ELSE
+ sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
+ ENDIF
+
+ IF ( transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN
+ sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
+ ELSE
+ sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
+ ENDIF
+
+ ENDIF
+
+ IF ( transports_3d(4,jsec,jseg,jk) .GE. 0.0 ) THEN
+ sec%transport(7,jclass) = sec%transport(7,jclass)+transports_3d(4,jsec,jseg,jk)
+ ELSE
+ sec%transport(8,jclass) = sec%transport(8,jclass)+transports_3d(4,jsec,jseg,jk)
+ ENDIF
+
+ IF ( transports_3d(5,jsec,jseg,jk) .GE. 0.0 ) THEN
+ sec%transport( 9,jclass) = sec%transport( 9,jclass)+transports_3d(5,jsec,jseg,jk)
+ ELSE
+ sec%transport(10,jclass) = sec%transport(10,jclass)+transports_3d(5,jsec,jseg,jk)
+ ENDIF
+
+ ELSE
+ sec%transport( 3,jclass) = 0._wp
+ sec%transport( 4,jclass) = 0._wp
+ sec%transport( 5,jclass) = 0._wp
+ sec%transport( 6,jclass) = 0._wp
+ sec%transport( 7,jclass) = 0._wp
+ sec%transport( 8,jclass) = 0._wp
+ sec%transport( 9,jclass) = 0._wp
+ sec%transport(10,jclass) = 0._wp
+ ENDIF
+
+ ENDIF ! end of test if point is in class
+
+ ENDDO ! end of loop on the classes
+
+ ENDDO ! loop over jk
+
+#if defined key_lim2 || defined key_lim3
+
+ !ICE CASE
+ IF( sec%ll_ice_section )THEN
+
+ IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN
+ sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg)
+ ELSE
+ sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg)
+ ENDIF
+
+ IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN
+ sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg)
+ ELSE
+ sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg)
+ ENDIF
+
+ ENDIF !end of ice case
+#endif
+
+ ENDDO !end of loop on the segment
+
+ ELSE !if sec%nb_point =0
+ sec%transport(1:2,:)=0.
+ IF (sec%llstrpond) sec%transport(3:10,:)=0.
+ IF (sec%ll_ice_section) sec%transport( 11:14,:)=0.
+ ENDIF !end of sec%nb_point =0 case
+
+ END SUBROUTINE dia_dct_sum
+
+ SUBROUTINE dia_dct_sum_h(sec,jsec)
+ !!-------------------------------------------------------------
+ !! Exactly as dia_dct_sum but for hourly files containing data summed at each time step
+ !!
+ !! Purpose: Average the transport over nn_dctwri time steps
+ !! and sum over the density/salinity/temperature/depth classes
+ !!
+ !! Method:
+ !! Sum over relevant grid cells to obtain values
+ !! for each
+ !! There are several loops:
+ !! loop on the segment between 2 nodes
+ !! loop on the level jk
+ !! loop on the density/temperature/salinity/level classes
+ !! test on the density/temperature/salinity/level
+ !!
+ !! ** Method :Transport through a given section is equal to the sum of transports
+ !! computed on each proc.
+ !! On each proc,transport is equal to the sum of transport computed through
+ !! segments linking each point of sec%listPoint with the next one.
+ !!
+ !!-------------------------------------------------------------
+ !! * arguments
+ TYPE(SECTION),INTENT(INOUT) :: sec
+ INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section
+
+ TYPE(POINT_SECTION) :: k
+ INTEGER :: jk,jseg,jclass !loop on level/segment/classes
+ REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
+ !!-------------------------------------------------------------
+
+ !! Sum the relevant segments to obtain values for each class
+ IF(sec%nb_point .NE. 0)THEN
+
+ !--------------------------------------!
+ ! LOOP ON THE SEGMENT BETWEEN 2 NODES !
+ !--------------------------------------!
+ DO jseg=1,MAX(sec%nb_point-1,0)
+
+ !-------------------------------------------------------------------------------------------
+ ! Select the appropriate coordinate for computing the velocity of the segment
+ !
+ ! CASE(0) Case (2)
+ ! ------- --------
+ ! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j)
+ ! F(i,j)----------V(i+1,j)-------F(i+1,j) |
+ ! |
+ ! |
+ ! |
+ ! Case (3) U(i,j)
+ ! -------- |
+ ! |
+ ! listPoint(jseg+1) F(i,j+1) |
+ ! | |
+ ! | |
+ ! | listPoint(jseg+1) F(i,j-1)
+ ! |
+ ! |
+ ! U(i,j+1)
+ ! | Case(1)
+ ! | ------
+ ! |
+ ! | listPoint(jseg+1) listPoint(jseg)
+ ! | F(i-1,j)-----------V(i,j) -------f(jseg)
+ ! listPoint(jseg) F(i,j)
+ !
+ !-------------------------------------------------------------------------------------------
+
+ SELECT CASE( sec%direction(jseg) )
+ CASE(0) ; k = sec%listPoint(jseg)
+ CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
+ CASE(2) ; k = sec%listPoint(jseg)
+ CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
+ END SELECT
+
+ !---------------------------|
+ ! LOOP ON THE LEVEL |
+ !---------------------------|
+ !Sum of the transport on the vertical
+ DO jk=1,mbathy(k%I,k%J)
+
+ ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
+ SELECT CASE( sec%direction(jseg) )
+ CASE(0,1)
+ ztn = interp(k%I,k%J,jk,'V',0)
+ zsn = interp(k%I,k%J,jk,'V',1)
+ zrhop = interp(k%I,k%J,jk,'V',2)
+ zrhoi = interp(k%I,k%J,jk,'V',3)
+ zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1)
+ CASE(2,3)
+ ztn = interp(k%I,k%J,jk,'U',0)
+ zsn = interp(k%I,k%J,jk,'U',1)
+ zrhop = interp(k%I,k%J,jk,'U',2)
+ zrhoi = interp(k%I,k%J,jk,'U',3)
zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1)
END SELECT
@@ -1472,15 +1712,19 @@
IF( sec%llstrpond )THEN
- IF ( transports_3d_h(2,jsec,jseg,jk) .GE. 0.0 ) THEN
- sec%transport_h(3,jclass) = sec%transport_h(3,jclass)+transports_3d_h(2,jsec,jseg,jk)
- ELSE
- sec%transport_h(4,jclass) = sec%transport_h(4,jclass)+transports_3d_h(2,jsec,jseg,jk)
- ENDIF
-
- IF ( transports_3d_h(3,jsec,jseg,jk) .GE. 0.0 ) THEN
- sec%transport_h(5,jclass) = sec%transport_h(5,jclass)+transports_3d_h(3,jsec,jseg,jk)
- ELSE
- sec%transport_h(6,jclass) = sec%transport_h(6,jclass)+transports_3d_h(3,jsec,jseg,jk)
- ENDIF
+ IF( transports_3d_h(1,jsec,jseg,jk) .NE. 0._wp ) THEN
+
+ IF (transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN
+ sec%transport_h(3,jclass) = sec%transport_h(3,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
+ ELSE
+ sec%transport_h(4,jclass) = sec%transport_h(4,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
+ ENDIF
+
+ IF ( transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN
+ sec%transport_h(5,jclass) = sec%transport_h(5,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
+ ELSE
+ sec%transport_h(6,jclass) = sec%transport_h(6,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
+ ENDIF
+
+ ENDIF
IF ( transports_3d_h(4,jsec,jseg,jk) .GE. 0.0 ) THEN
@@ -1572,12 +1816,36 @@
!
REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
+ CHARACTER(len=3) :: noos_sect_name ! Classname
+ CHARACTER(len=25) :: noos_var_sect_name
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: noos_iom_dummy
+ INTEGER :: IERR
+
+ REAL(wp), DIMENSION(3) :: tmp_iom_output
+ REAL(wp) :: max_iom_val
+
!!-------------------------------------------------------------
+
+
+
+ IF( lwp ) THEN
+ WRITE(numout,*) " "
+ WRITE(numout,*) "dia_dct_wri_NOOS: write transports through sections at timestep: ", kt
+ WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
+ ENDIF
+
CALL wrk_alloc(nb_type_class , zsumclasses )
zsumclasses(:)=0._wp
zslope = sec%slopeSection
-
- WRITE(numdct_NOOS,'(I4,a1,I2,a1,I2,a12,i3,a17,i3,a10,a25)') nyear,'.',nmonth,'.',nday,' Transect:',ksec-1,' No. of layers:',sec%nb_class-1,' Name: ',sec%name
-
+
+ IF( lwp ) THEN
+ IF ( ln_dct_ascii ) THEN
+ WRITE(numdct_NOOS,'(I4,a1,I2,a1,I2,a12,i3,a17,i3,a10,a25)') nyear,'.',nmonth,'.',nday,' Transect:',ksec-1,' No. of layers:',sec%nb_class-1,' Name: ',sec%name
+ ELSE
+ WRITE(numdct_NOOS) nyear,nmonth,nday,ksec-1,sec%nb_class-1,sec%name
+ ENDIF
+ ENDIF
+
+ ! Sum all classes together, to give one values per type (pos tran, neg vol trans etc...).
DO jclass=1,MAX(1,sec%nb_class-1)
zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
@@ -1587,14 +1855,407 @@
zbnd1 = 0._wp
zbnd2 = 0._wp
-
- IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
- WRITE(numdct_NOOS,'(9e12.4E2)') -(zsumclasses( 1)+zsumclasses( 2)), -zsumclasses( 2),-zsumclasses( 1), &
- -(zsumclasses( 7)+zsumclasses( 8)), -zsumclasses( 8),-zsumclasses( 7), &
- -(zsumclasses( 9)+zsumclasses(10)), -zsumclasses(10),-zsumclasses( 9)
- ELSE
- WRITE(numdct_NOOS,'(9e12.4E2)') zsumclasses( 1)+zsumclasses( 2) , zsumclasses( 1), zsumclasses( 2), &
- zsumclasses( 7)+zsumclasses( 8) , zsumclasses( 7), zsumclasses( 8), &
- zsumclasses( 9)+zsumclasses(10) , zsumclasses( 9), zsumclasses(10)
- ENDIF
+
+
+
+ write (noos_sect_name, "(I0.3)") ksec
+
+ IF ( ln_dct_iom_cont ) THEN
+ max_iom_val = 1.e10
+ ALLOCATE( noos_iom_dummy(jpi,jpj,3), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_dct_wri_NOOS: failed to allocate noos_iom_dummy array' )
+ ENDIF
+
+! JT I think changing the sign on the output based on the zslope value is redunant.
+!
+! IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
+!
+! IF( lwp ) THEN
+! WRITE(numdct_NOOS,'(9e12.4E2)') -(zsumclasses( 1)+zsumclasses( 2)), -zsumclasses( 2),-zsumclasses( 1), &
+! -(zsumclasses( 7)+zsumclasses( 8)), -zsumclasses( 8),-zsumclasses( 7), &
+! -(zsumclasses( 9)+zsumclasses(10)), -zsumclasses(10),-zsumclasses( 9)
+! CALL FLUSH(numdct_NOOS)
+! endif
+
+!
+! IF ( ln_dct_iom_cont ) THEN
+!
+! noos_iom_dummy(:,:,:) = 0.
+!
+! tmp_iom_output(:) = 0.
+! tmp_iom_output(1) = -(zsumclasses( 1)+zsumclasses( 2))
+! tmp_iom_output(2) = -zsumclasses( 2)
+! tmp_iom_output(3) = -zsumclasses( 1)
+!
+! ! Convert to Sv
+! tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
+! tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
+! tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
+!
+! ! limit maximum and minimum values in iom_put
+! if ( tmp_iom_output(1) .gt. max_iom_val ) tmp_iom_output(1) = max_iom_val
+! if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+! if ( tmp_iom_output(2) .gt. max_iom_val ) tmp_iom_output(2) = max_iom_val
+! if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+! if ( tmp_iom_output(3) .gt. max_iom_val ) tmp_iom_output(3) = max_iom_val
+! if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+!
+! ! Set NaN's to Zero
+! if ( tmp_iom_output(1) .ne. tmp_iom_output(1) ) tmp_iom_output(1) = max_iom_val*2
+! if ( tmp_iom_output(2) .ne. tmp_iom_output(2) ) tmp_iom_output(1) = max_iom_val*2
+! if ( tmp_iom_output(3) .ne. tmp_iom_output(3) ) tmp_iom_output(1) = max_iom_val*2
+!
+! noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+! noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+! noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+!
+! noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
+! if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+! CALL iom_put( noos_var_sect_name, noos_iom_dummy )
+! noos_iom_dummy(:,:,:) = 0.
+! tmp_iom_output(:) = 0.
+! tmp_iom_output(1) = -(zsumclasses( 7)+zsumclasses( 8))
+! tmp_iom_output(2) = -zsumclasses( 8)
+! tmp_iom_output(3) = -zsumclasses( 7)
+!
+! ! Convert to TJ/s
+! tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
+! tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
+! tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
+!
+! ! limit maximum and minimum values in iom_put
+! if ( tmp_iom_output(1) .gt. max_iom_val ) tmp_iom_output(1) = max_iom_val
+! if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+! if ( tmp_iom_output(2) .gt. max_iom_val ) tmp_iom_output(2) = max_iom_val
+! if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+! if ( tmp_iom_output(3) .gt. max_iom_val ) tmp_iom_output(3) = max_iom_val
+! if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+!
+! ! Set NaN's to Zero
+! if ( tmp_iom_output(1) .ne. tmp_iom_output(1) ) tmp_iom_output(1) = max_iom_val*2
+! if ( tmp_iom_output(2) .ne. tmp_iom_output(2) ) tmp_iom_output(1) = max_iom_val*2
+! if ( tmp_iom_output(3) .ne. tmp_iom_output(3) ) tmp_iom_output(1) = max_iom_val*2
+!
+! noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+! noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+! noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+!
+! !noos_iom_dummy(:,:,1) = -(zsumclasses( 7)+zsumclasses( 8))
+! !noos_iom_dummy(:,:,2) = -zsumclasses( 8)
+! !noos_iom_dummy(:,:,3) = -zsumclasses( 7)
+!
+! noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_heat'
+! if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+! CALL iom_put( noos_var_sect_name, noos_iom_dummy )
+!
+! noos_iom_dummy(:,:,:) = 0.
+! tmp_iom_output(:) = 0.
+! tmp_iom_output(1) = -(zsumclasses( 9)+zsumclasses( 10))
+! tmp_iom_output(2) = -zsumclasses( 10)
+! tmp_iom_output(3) = -zsumclasses( 9)
+!
+! ! Convert to MT/s
+! tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
+! tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
+! tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
+!
+! ! limit maximum and minimum values in iom_put
+! if ( tmp_iom_output(1) .gt. max_iom_val ) tmp_iom_output(1) = max_iom_val
+! if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+! if ( tmp_iom_output(2) .gt. max_iom_val ) tmp_iom_output(2) = max_iom_val
+! if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+! if ( tmp_iom_output(3) .gt. max_iom_val ) tmp_iom_output(3) = max_iom_val
+! if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+!
+! ! Set NaN's to Zero
+! if ( tmp_iom_output(1) .ne. tmp_iom_output(1) ) tmp_iom_output(1) = max_iom_val*2
+! if ( tmp_iom_output(2) .ne. tmp_iom_output(2) ) tmp_iom_output(1) = max_iom_val*2
+! if ( tmp_iom_output(3) .ne. tmp_iom_output(3) ) tmp_iom_output(1) = max_iom_val*2
+!
+! noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+! noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+! noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+!
+! !noos_iom_dummy(:,:,1) = -(zsumclasses( 9)+zsumclasses( 10))
+! !noos_iom_dummy(:,:,2) = -zsumclasses( 10)
+! !noos_iom_dummy(:,:,3) = -zsumclasses( 9)
+!
+! noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_salt'
+! if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+! CALL iom_put( noos_var_sect_name, noos_iom_dummy )
+! noos_iom_dummy(:,:,:) = 0.
+! tmp_iom_output(:) = 0.
+! ENDIF
+! ELSE
+! IF( lwp ) THEN
+! WRITE(numdct_NOOS,'(9e12.4E2)') zsumclasses( 1)+zsumclasses( 2) , zsumclasses( 1), zsumclasses( 2), &
+! zsumclasses( 7)+zsumclasses( 8) , zsumclasses( 7), zsumclasses( 8), &
+! zsumclasses( 9)+zsumclasses(10) , zsumclasses( 9), zsumclasses(10)
+! CALL FLUSH(numdct_NOOS)
+! endif
+!
+!
+! IF ( ln_dct_iom_cont ) THEN
+!
+! noos_iom_dummy(:,:,:) = 0.
+! tmp_iom_output(:) = 0.
+!
+! tmp_iom_output(1) = (zsumclasses( 1)+zsumclasses( 2))
+! tmp_iom_output(2) = zsumclasses( 1)
+! tmp_iom_output(3) = zsumclasses( 2)
+!
+! ! Convert to Sv
+! tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
+! tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
+! tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
+!
+! ! limit maximum and minimum values in iom_put
+! if ( tmp_iom_output(1) .gt. max_iom_val ) tmp_iom_output(1) = max_iom_val
+! if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+! if ( tmp_iom_output(2) .gt. max_iom_val ) tmp_iom_output(2) = max_iom_val
+! if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+! if ( tmp_iom_output(3) .gt. max_iom_val ) tmp_iom_output(3) = max_iom_val
+! if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+!
+! ! Set NaN's to Zero
+! if ( tmp_iom_output(1) .ne. tmp_iom_output(1) ) tmp_iom_output(1) = max_iom_val*2
+! if ( tmp_iom_output(2) .ne. tmp_iom_output(2) ) tmp_iom_output(1) = max_iom_val*2
+! if ( tmp_iom_output(3) .ne. tmp_iom_output(3) ) tmp_iom_output(1) = max_iom_val*2
+!
+! noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+! noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+! noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+!
+! !noos_iom_dummy(:,:,1) = (zsumclasses( 1)+zsumclasses( 2))
+! !noos_iom_dummy(:,:,2) = zsumclasses( 1)
+! !noos_iom_dummy(:,:,3) = zsumclasses( 2)
+!
+!
+!
+! noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
+! if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+! CALL iom_put( noos_var_sect_name, noos_iom_dummy )
+! noos_iom_dummy(:,:,:) = 0.
+! tmp_iom_output(:) = 0.
+!
+! tmp_iom_output(1) = (zsumclasses( 7)+zsumclasses( 8))
+! tmp_iom_output(2) = zsumclasses( 7)
+! tmp_iom_output(3) = zsumclasses( 8)
+!
+! ! Convert to TJ/s
+! tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
+! tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
+! tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
+!
+! ! limit maximum and minimum values in iom_put
+! if ( tmp_iom_output(1) .gt. max_iom_val ) tmp_iom_output(1) = max_iom_val
+! if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+! if ( tmp_iom_output(2) .gt. max_iom_val ) tmp_iom_output(2) = max_iom_val
+! if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+! if ( tmp_iom_output(3) .gt. max_iom_val ) tmp_iom_output(3) = max_iom_val
+! if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+!
+! ! Set NaN's to Zero
+! if ( tmp_iom_output(1) .ne. tmp_iom_output(1) ) tmp_iom_output(1) = max_iom_val*2
+! if ( tmp_iom_output(2) .ne. tmp_iom_output(2) ) tmp_iom_output(1) = max_iom_val*2
+! if ( tmp_iom_output(3) .ne. tmp_iom_output(3) ) tmp_iom_output(1) = max_iom_val*2
+!
+! noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+! noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+! noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+!
+! !noos_iom_dummy(:,:,1) = (zsumclasses( 7)+zsumclasses( 8))
+! !noos_iom_dummy(:,:,2) = zsumclasses( 7)
+! !noos_iom_dummy(:,:,3) = zsumclasses( 8)
+!
+! noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_heat'
+! if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+! CALL iom_put(noos_var_sect_name, noos_iom_dummy )
+! noos_iom_dummy(:,:,:) = 0.
+! tmp_iom_output(:) = 0.
+!
+! tmp_iom_output(1) = (zsumclasses( 9)+zsumclasses( 10))
+! tmp_iom_output(2) = zsumclasses( 9)
+! tmp_iom_output(3) = zsumclasses( 10)
+!
+! ! Convert to MT/s
+! tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
+! tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
+! tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
+!
+!
+! ! limit maximum and minimum values in iom_put
+! if ( tmp_iom_output(1) .gt. max_iom_val ) tmp_iom_output(1) = max_iom_val
+! if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+! if ( tmp_iom_output(2) .gt. max_iom_val ) tmp_iom_output(2) = max_iom_val
+! if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+! if ( tmp_iom_output(3) .gt. max_iom_val ) tmp_iom_output(3) = max_iom_val
+! if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+!
+! ! Set NaN's to Zero
+! if ( tmp_iom_output(1) .ne. tmp_iom_output(1) ) tmp_iom_output(1) = max_iom_val*2
+! if ( tmp_iom_output(2) .ne. tmp_iom_output(2) ) tmp_iom_output(1) = max_iom_val*2
+! if ( tmp_iom_output(3) .ne. tmp_iom_output(3) ) tmp_iom_output(1) = max_iom_val*2
+!
+! noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+! noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+! noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+!
+! !noos_iom_dummy(:,:,1) = (zsumclasses( 9)+zsumclasses( 10))
+! !noos_iom_dummy(:,:,2) = zsumclasses( 9)
+! !noos_iom_dummy(:,:,3) = zsumclasses( 10)
+!
+! noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_salt'
+! if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+! CALL iom_put(noos_var_sect_name, noos_iom_dummy )
+! noos_iom_dummy(:,:,:) = 0.
+! tmp_iom_output(:) = 0.
+! ENDIF
+!
+! ENDIF
+
+
+
+
+
+
+
+
+
+
+ IF( lwp ) THEN
+ IF ( ln_dct_ascii ) THEN
+ !WRITE(numdct_NOOS,'(9e12.4E2)') zsumclasses( 1)+zsumclasses( 2) , zsumclasses( 1), zsumclasses( 2), &
+ WRITE(numdct_NOOS,'(3F18.3,6e16.8E2)') zsumclasses( 1)+zsumclasses( 2) , zsumclasses( 1), zsumclasses( 2), &
+ zsumclasses( 7)+zsumclasses( 8) , zsumclasses( 7), zsumclasses( 8), &
+ zsumclasses( 9)+zsumclasses(10) , zsumclasses( 9), zsumclasses(10)
+ CALL FLUSH(numdct_NOOS)
+ ELSE
+ WRITE(numdct_NOOS) zsumclasses( 1)+zsumclasses( 2) , zsumclasses( 1), zsumclasses( 2), &
+ zsumclasses( 7)+zsumclasses( 8) , zsumclasses( 7), zsumclasses( 8), &
+ zsumclasses( 9)+zsumclasses(10) , zsumclasses( 9), zsumclasses(10)
+ CALL FLUSH(numdct_NOOS)
+ ENDIF
+ ENDIF
+
+ IF ( ln_dct_iom_cont ) THEN
+
+ noos_iom_dummy(:,:,:) = 0.
+ tmp_iom_output(:) = 0.
+
+ tmp_iom_output(1) = (zsumclasses( 1)+zsumclasses( 2))
+ tmp_iom_output(2) = zsumclasses( 1)
+ tmp_iom_output(3) = zsumclasses( 2)
+
+ ! Convert to Sv
+ tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
+ tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
+ tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
+
+ ! limit maximum and minimum values in iom_put
+ if ( tmp_iom_output(1) .gt. max_iom_val ) tmp_iom_output(1) = max_iom_val
+ if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+ if ( tmp_iom_output(2) .gt. max_iom_val ) tmp_iom_output(2) = max_iom_val
+ if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+ if ( tmp_iom_output(3) .gt. max_iom_val ) tmp_iom_output(3) = max_iom_val
+ if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+
+ ! Set NaN's to Zero
+ if ( tmp_iom_output(1) .ne. tmp_iom_output(1) ) tmp_iom_output(1) = max_iom_val*2
+ if ( tmp_iom_output(2) .ne. tmp_iom_output(2) ) tmp_iom_output(1) = max_iom_val*2
+ if ( tmp_iom_output(3) .ne. tmp_iom_output(3) ) tmp_iom_output(1) = max_iom_val*2
+
+ noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+ noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+ noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+
+ !noos_iom_dummy(:,:,1) = (zsumclasses( 1)+zsumclasses( 2))
+ !noos_iom_dummy(:,:,2) = zsumclasses( 1)
+ !noos_iom_dummy(:,:,3) = zsumclasses( 2)
+
+
+
+ noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
+ if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+ CALL iom_put( noos_var_sect_name, noos_iom_dummy )
+ noos_iom_dummy(:,:,:) = 0.
+ tmp_iom_output(:) = 0.
+
+ tmp_iom_output(1) = (zsumclasses( 7)+zsumclasses( 8))
+ tmp_iom_output(2) = zsumclasses( 7)
+ tmp_iom_output(3) = zsumclasses( 8)
+
+ ! Convert to TJ/s
+ tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
+ tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
+ tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
+
+ ! limit maximum and minimum values in iom_put
+ if ( tmp_iom_output(1) .gt. max_iom_val ) tmp_iom_output(1) = max_iom_val
+ if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+ if ( tmp_iom_output(2) .gt. max_iom_val ) tmp_iom_output(2) = max_iom_val
+ if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+ if ( tmp_iom_output(3) .gt. max_iom_val ) tmp_iom_output(3) = max_iom_val
+ if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+
+ ! Set NaN's to Zero
+ if ( tmp_iom_output(1) .ne. tmp_iom_output(1) ) tmp_iom_output(1) = max_iom_val*2
+ if ( tmp_iom_output(2) .ne. tmp_iom_output(2) ) tmp_iom_output(1) = max_iom_val*2
+ if ( tmp_iom_output(3) .ne. tmp_iom_output(3) ) tmp_iom_output(1) = max_iom_val*2
+
+ noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+ noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+ noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+
+ !noos_iom_dummy(:,:,1) = (zsumclasses( 7)+zsumclasses( 8))
+ !noos_iom_dummy(:,:,2) = zsumclasses( 7)
+ !noos_iom_dummy(:,:,3) = zsumclasses( 8)
+
+ noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_heat'
+ if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+ CALL iom_put(noos_var_sect_name, noos_iom_dummy )
+ noos_iom_dummy(:,:,:) = 0.
+ tmp_iom_output(:) = 0.
+
+ tmp_iom_output(1) = (zsumclasses( 9)+zsumclasses( 10))
+ tmp_iom_output(2) = zsumclasses( 9)
+ tmp_iom_output(3) = zsumclasses( 10)
+
+ ! Convert to MT/s
+ tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
+ tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
+ tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
+
+
+ ! limit maximum and minimum values in iom_put
+ if ( tmp_iom_output(1) .gt. max_iom_val ) tmp_iom_output(1) = max_iom_val
+ if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+ if ( tmp_iom_output(2) .gt. max_iom_val ) tmp_iom_output(2) = max_iom_val
+ if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+ if ( tmp_iom_output(3) .gt. max_iom_val ) tmp_iom_output(3) = max_iom_val
+ if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+
+ ! Set NaN's to Zero
+ if ( tmp_iom_output(1) .ne. tmp_iom_output(1) ) tmp_iom_output(1) = max_iom_val*2
+ if ( tmp_iom_output(2) .ne. tmp_iom_output(2) ) tmp_iom_output(1) = max_iom_val*2
+ if ( tmp_iom_output(3) .ne. tmp_iom_output(3) ) tmp_iom_output(1) = max_iom_val*2
+
+ noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+ noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+ noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+
+ !noos_iom_dummy(:,:,1) = (zsumclasses( 9)+zsumclasses( 10))
+ !noos_iom_dummy(:,:,2) = zsumclasses( 9)
+ !noos_iom_dummy(:,:,3) = zsumclasses( 10)
+
+ noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_salt'
+ if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+ CALL iom_put(noos_var_sect_name, noos_iom_dummy )
+ noos_iom_dummy(:,:,:) = 0.
+ tmp_iom_output(:) = 0.
+
+
+ DEALLOCATE(noos_iom_dummy)
+ ENDIF
+
DO jclass=1,MAX(1,sec%nb_class-1)
@@ -1641,15 +2302,29 @@
!write volume transport per class
- IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
- WRITE(numdct_NOOS,'(9e12.4E2)') -(sec%transport( 1,jclass)+sec%transport( 2,jclass)),-sec%transport( 2,jclass),-sec%transport( 1,jclass), &
- -(sec%transport( 7,jclass)+sec%transport( 8,jclass)),-sec%transport( 8,jclass),-sec%transport( 7,jclass), &
- -(sec%transport( 9,jclass)+sec%transport(10,jclass)),-sec%transport(10,jclass),-sec%transport( 9,jclass)
- ELSE
- WRITE(numdct_NOOS,'(9e12.4E2)') sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
- sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
- sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
+ IF( lwp ) THEN
+
+ IF ( ln_dct_ascii ) THEN
+ CALL FLUSH(numdct_NOOS)
+
+ !WRITE(numdct_NOOS,'(9e12.4E2)') sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
+ ! sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
+ ! sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
+ WRITE(numdct_NOOS,'(3F18.3,6e16.8E2)') sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
+ sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
+ sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
+ ELSE
+
+ CALL FLUSH(numdct_NOOS)
+ WRITE(numdct_NOOS) sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
+ sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
+ sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
+ ENDIF
ENDIF
ENDDO
+
+ !IF ( ln_dct_ascii ) THEN
+ if ( lwp ) CALL FLUSH(numdct_NOOS)
+ !ENDIF
CALL wrk_dealloc(nb_type_class , zsumclasses )
@@ -1671,5 +2346,5 @@
!!-------------------------------------------------------------
!!arguments
- INTEGER, INTENT(IN) :: hr ! hour
+ INTEGER, INTENT(IN) :: hr ! hour => effectively kt/12
TYPE(SECTION), INTENT(INOUT) :: sec ! section to write
INTEGER ,INTENT(IN) :: ksec ! section number
@@ -1682,38 +2357,88 @@
!
REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
+ CHARACTER(len=3) :: noos_sect_name ! Classname
+ CHARACTER(len=25) :: noos_var_sect_name
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: noos_iom_dummy
+ INTEGER :: IERR
+
!!-------------------------------------------------------------
+ IF( lwp ) THEN
+ WRITE(numout,*) " "
+ WRITE(numout,*) "dia_dct_wri_NOOS_h: write transports through section Transect:",ksec-1," at timestep: ", hr
+ WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
+ ENDIF
+
CALL wrk_alloc(nb_type_class , zsumclasses )
+
+
+ write (noos_sect_name, "(I03)") ksec
+
+ ALLOCATE( noos_iom_dummy(jpi,jpj,3), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_dct_wri_NOOS_h: failed to allocate noos_iom_dummy array' )
+
+
+
+
zsumclasses(:)=0._wp
zslope = sec%slopeSection
+ ! Sum up all classes, to give the total per type (pos vol trans, neg vol trans etc...)
DO jclass=1,MAX(1,sec%nb_class-1)
zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport_h(1:nb_type_class,jclass)
ENDDO
-
+
+
+ ! JT I think changing the sign of output according to the zslope is redundant
+
!write volume transport per class
- IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
- z_hr_output(ksec,hr,1)=-(zsumclasses(1)+zsumclasses(2))
- ELSE
- z_hr_output(ksec,hr,1)= (zsumclasses(1)+zsumclasses(2))
- ENDIF
-
+ ! Sum positive and vol trans for all classes in first cell of array
+
+ z_hr_output(ksec,1,1)= (zsumclasses(1)+zsumclasses(2))
+ z_hr_output(ksec,2,1)= zsumclasses(1)
+ z_hr_output(ksec,3,1)= zsumclasses(2)
+
+ ! Sum positive and vol trans for each classes in following cell of array
DO jclass=1,MAX(1,sec%nb_class-1)
- IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
- z_hr_output(ksec,hr,jclass+1)=-(sec%transport_h(1,jclass)+sec%transport_h(2,jclass))
+ z_hr_output(ksec,1,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass))
+ z_hr_output(ksec,2,jclass+1)= sec%transport_h(1,jclass)
+ z_hr_output(ksec,3,jclass+1)= sec%transport_h(2,jclass)
+ ENDDO
+
+
+ IF( lwp ) THEN
+ ! JT IF ( hr .eq. 48._wp ) THEN
+ ! JT WRITE(numdct_NOOS_h,'(I4,a1,I2,a1,I2,a12,i3,a17,i3)') nyear,'.',nmonth,'.',nday,' Transect:',ksec-1,' No. of layers:',sec%nb_class-1
+ ! JT DO jhr=25,48
+ ! JT WRITE(numdct_NOOS_h,'(11F12.1)') z_hr_output(ksec,jhr,1), (z_hr_output(ksec,jhr,jclass+1), jclass=1,MAX(1,10) )
+ ! JT ENDDO
+ ! JT ENDIF
+
+
+
+ IF ( ln_dct_ascii ) THEN
+ WRITE(numdct_NOOS_h,'(I4,a1,I2,a1,I2,a1,I2,a1,I2,a12,i3,a17,i3)') nyear,'.',nmonth,'.',nday,'.',MOD(hr,24),'.',0,' Transect:',ksec-1,' No. of layers:',sec%nb_class-1
+ WRITE(numdct_NOOS_h,'(11F18.3)') z_hr_output(ksec,1,1), (z_hr_output(ksec,1,jclass+1), jclass=1,MAX(1,10) )
+ WRITE(numdct_NOOS_h,'(11F18.3)') z_hr_output(ksec,2,1), (z_hr_output(ksec,2,jclass+1), jclass=1,MAX(1,10) )
+ WRITE(numdct_NOOS_h,'(11F18.3)') z_hr_output(ksec,3,1), (z_hr_output(ksec,3,jclass+1), jclass=1,MAX(1,10) )
+ CALL FLUSH(numdct_NOOS_h)
ELSE
- z_hr_output(ksec,hr,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass))
+ WRITE(numdct_NOOS_h) nyear,nmonth,nday,MOD(hr,24),ksec-1,sec%nb_class-1
+ WRITE(numdct_NOOS_h) z_hr_output(ksec,1,1), (z_hr_output(ksec,1,jclass+1), jclass=1,MAX(1,10) )
+ WRITE(numdct_NOOS_h) z_hr_output(ksec,2,1), (z_hr_output(ksec,2,jclass+1), jclass=1,MAX(1,10) )
+ WRITE(numdct_NOOS_h) z_hr_output(ksec,3,1), (z_hr_output(ksec,3,jclass+1), jclass=1,MAX(1,10) )
+ CALL FLUSH(numdct_NOOS_h)
ENDIF
- ENDDO
-
- IF ( hr .eq. 48._wp ) THEN
- WRITE(numdct_NOOS_h,'(I4,a1,I2,a1,I2,a12,i3,a17,i3)') nyear,'.',nmonth,'.',nday,' Transect:',ksec-1,' No. of layers:',sec%nb_class-1
- DO jhr=25,48
- WRITE(numdct_NOOS_h,'(11F12.1)') z_hr_output(ksec,jhr,1), (z_hr_output(ksec,jhr,jclass+1), jclass=1,MAX(1,10) )
- ENDDO
- ENDIF
+
+
+ ENDIF
+
CALL wrk_dealloc(nb_type_class , zsumclasses )
+
+ DEALLOCATE(noos_iom_dummy)
+
+
END SUBROUTINE dia_dct_wri_NOOS_h
@@ -1730,8 +2455,8 @@
!!
!! 2. Write heat transports in "heat_transport"
- !! Unit: Peta W : area * Velocity * T * rhop * Cp * 1.e-15
+ !! Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15
!!
!! 3. Write salt transports in "salt_transport"
- !! Unit: 10^9 Kg/m^2/s : area * Velocity * S * rhop * 1.e-9
+ !! Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6
!!
!!-------------------------------------------------------------
@@ -1810,11 +2535,11 @@
WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
jclass,classe,zbnd1,zbnd2,&
- sec%transport(7,jclass)*1.e-15,sec%transport(8,jclass)*1.e-15, &
- ( sec%transport(7,jclass)+sec%transport(8,jclass) )*1.e-15
+ sec%transport(7,jclass)*1000._wp*rcp/1.e15,sec%transport(8,jclass)*1000._wp*rcp/1.e15, &
+ ( sec%transport(7,jclass)+sec%transport(8,jclass) )*1000._wp*rcp/1.e15
!write salt transport per class
WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
jclass,classe,zbnd1,zbnd2,&
- sec%transport(9,jclass)*1.e-9,sec%transport(10,jclass)*1.e-9,&
- (sec%transport(9,jclass)+sec%transport(10,jclass))*1.e-9
+ sec%transport(9,jclass)*1000._wp/1.e9,sec%transport(10,jclass)*1000._wp/1.e9,&
+ (sec%transport(9,jclass)+sec%transport(10,jclass))*1000._wp/1.e9
ENDIF
@@ -1835,11 +2560,11 @@
WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
jclass,"total",zbnd1,zbnd2,&
- zsumclasses(7)*1.e-15,zsumclasses(8)*1.e-15,&
- (zsumclasses(7)+zsumclasses(8) )*1.e-15
+ zsumclasses(7)* 1000._wp*rcp/1.e15,zsumclasses(8)* 1000._wp*rcp/1.e15,&
+ (zsumclasses(7)+zsumclasses(8) )* 1000._wp*rcp/1.e15
!write total salt transport
WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
jclass,"total",zbnd1,zbnd2,&
- zsumclasses(9)*1.e-9,zsumclasses(10)*1.e-9,&
- (zsumclasses(9)+zsumclasses(10))*1.e-9
+ zsumclasses(9)*1000._wp/1.e9,zsumclasses(10)*1000._wp/1.e9,&
+ (zsumclasses(9)+zsumclasses(10))*1000._wp/1.e9
ENDIF
@@ -1878,5 +2603,5 @@
!! | I | I+1 | Z=temperature/salinity/density at U-poinT
!! | | |
- !! ---------------------------------------- 1. Veritcal interpolation: compute zbis
+ !! ---------------------------------------- 1. Veritcale interpolation: compute zbis
!! | | | interpolation between ptab(I,J,K) and ptab(I,J,K+1)
!! | | | zbis =
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diapea.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diapea.F90 (revision 10958)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diapea.F90 (revision 10958)
@@ -0,0 +1,299 @@
+MODULE diapea
+ !!======================================================================
+ !! *** MODULE diapea ***
+ !! Potential Energy Anomaly
+ !!======================================================================
+ !! History : 3.6 ! 12/2016 (J Tinker) Original code
+ !!----------------------------------------------------------------------
+ USE oce ! ocean dynamics and tracers variables
+ USE dom_oce ! ocean space and time domain
+ USE in_out_manager ! I/O units
+ USE iom ! I/0 library
+ USE wrk_nemo ! working arrays
+ USE eosbn2 ! Equation of state - in situ and potential density
+ USE phycst ! physical constant
+
+ IMPLICIT NONE
+ PRIVATE
+
+ PUBLIC dia_pea_init ! routine called by nemogcm.F90
+ PUBLIC dia_pea ! routine called by diawri.F90
+ REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: pea,peat,peas
+ REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: wgt_co_mat ! Weighting array for proportion of grid shallower than cut off depth
+ REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: t_zmean, s_zmean !Depth mean temperature and salinity: 2d fields
+ REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: t_zmean_mat, s_zmean_mat !Depth mean temperature and salinity: 3d fields
+ REAL(wp) :: zcutoff
+ LOGICAL , PUBLIC :: ln_pea ! region mean calculation
+ !!----------------------------------------------------------------------
+ !! NEMO/OPA 3.6 , NEMO Consortium (2014)
+ !! $Id$
+ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+ !!----------------------------------------------------------------------
+CONTAINS
+
+ SUBROUTINE dia_pea_init
+ ! Local variables
+ INTEGER :: ji,jj,jk ! Dummy loop indices
+ REAL(wp) :: sumz,tmpsumz
+ INTEGER :: ierr ! error integer for IOM_get
+ INTEGER :: ios ! Local integer output status for namelist read
+
+ zcutoff = 200.!200m
+
+ NAMELIST/nam_pea/ ln_pea
+
+
+ !
+ REWIND ( numnam_ref ) ! Read Namelist nam_diatmb in referdiatmbence namelist : TMB diagnostics
+ READ ( numnam_ref, nam_pea, IOSTAT=ios, ERR= 901 )
+901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_pea in reference namelist', lwp )
+
+ REWIND( numnam_cfg ) ! Namelist nam_diatmb in configuration namelist TMB diagnostics
+ READ ( numnam_cfg, nam_pea, IOSTAT = ios, ERR = 902 )
+902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_pea in configuration namelist', lwp )
+ IF(lwm) WRITE ( numond, nam_pea )
+
+ IF(lwp) THEN ! Control print
+ WRITE(numout,*)
+ WRITE(numout,*) 'dia_pea_init : Output potential energy anomaly Diagnostics'
+ WRITE(numout,*) '~~~~~~~~~~~~'
+ WRITE(numout,*) 'Namelist nam_pea : set pea output '
+ WRITE(numout,*) 'Switch for pea diagnostics (T) or not (F) ln_diaregmean = ', ln_pea
+ ENDIF
+
+
+ ALLOCATE( pea(jpi,jpj), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate pea array' )
+
+ ALLOCATE( peat(jpi,jpj), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate peat array' )
+
+ ALLOCATE( peas(jpi,jpj), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate peas array' )
+
+ ALLOCATE( t_zmean_mat(jpi,jpj,jpk), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate t_zmean_mat array' )
+
+ ALLOCATE( s_zmean_mat(jpi,jpj,jpk), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate s_zmean_mat array' )
+
+ ALLOCATE( t_zmean(jpi,jpj), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate t_zmean array' )
+
+ ALLOCATE( s_zmean(jpi,jpj), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate s_zmean array' )
+
+ ALLOCATE( wgt_co_mat(jpi,jpj,jpk), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate wgt_co_mat array' )
+
+
+ pea(:,:) = 0.
+ peat(:,:) = 0.
+ peas(:,:) = 0.
+
+ if ( ln_pea ) THEN
+
+ ! create wgt_co_mat mat, with the proportion of the grid (gdept_0) below cut off (200m)
+
+ DO jj = 1,jpj
+ DO ji = 1,jpi
+ IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
+ sumz = 0.
+ DO jk = 1,jpk
+ IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN
+ tmpsumz = sumz + e3t_n(ji,jj,jk)
+ IF (sumz .ge. zcutoff) THEN
+ ! Already too deep
+ wgt_co_mat(ji,jj,jk) = 0.
+ ELSE IF (tmpsumz .le. zcutoff) THEN
+ ! Too shallow
+ wgt_co_mat(ji,jj,jk) = 1.
+ ELSE
+ !proprotion of grid box above cut off depth
+ wgt_co_mat(ji,jj,jk) = (zcutoff-Sumz)/e3t_n(ji,jj,jk)
+ END IF
+ sumz = tmpsumz
+ endif
+ END DO
+
+ ELSE
+ !if land, set to 0.
+ DO jk = 1,jpk
+ wgt_co_mat(ji,jj,jk) = 0.
+ END DO
+
+ ENDIF
+ END DO
+ END DO
+ ENDIF
+
+
+ END SUBROUTINE dia_pea_init
+
+
+ SUBROUTINE dia_pea(kt)
+ INTEGER, INTENT( in ) :: kt ! ocean time-step index
+
+ INTEGER :: ji,jj,jk,ierr ! Dummy loop indices
+ REAL(wp) :: tmpdenom, tmpnum, maxz
+ !rau0
+
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ts_pea_mat,ts_pea_mat_TS_mean,ts_pea_mat_S_mean
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: tmp_pea_rho,tmp_pea_TS_mean_rho,tmp_pea_S_mean_rho
+ REAL(wp) :: int_y_pea,int_y_pea_t
+
+
+
+ ALLOCATE( ts_pea_mat(jpi,jpj,jpk,2), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate ts_pea_mat array' )
+ ALLOCATE( ts_pea_mat_TS_mean(jpi,jpj,jpk,2), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate ts_pea_mat_TS_mean array' )
+ ALLOCATE( ts_pea_mat_S_mean(jpi,jpj,jpk,2), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate ts_pea_mat_S_mean array' )
+ ALLOCATE( tmp_pea_rho(jpi,jpj,jpk), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate tmp_pea_rho array' )
+ ALLOCATE( tmp_pea_TS_mean_rho(jpi,jpj,jpk), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate tmp_pea_TS_mean_rho array' )
+ ALLOCATE( tmp_pea_S_mean_rho(jpi,jpj,jpk), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_pea_init: failed to allocate tmp_pea_S_mean_rho array' )
+
+ pea(:,:)=0.0d0
+ peaT(:,:)=0.0d0
+ peaS(:,:)=0.0d0
+
+
+
+ !calculate the depth mean temperature and salinity of the upper 200m. Save this into a 3d array. Set the value where tmask=0 to be tsn.
+
+ DO jj = 1,jpj
+ DO ji = 1,jpi
+ IF ( tmask(ji,jj,1) == 1.0_wp ) THEN ! if a sea point.
+
+ !Depth mean temperature
+ tmpdenom = 0.
+ tmpnum = 0.
+ DO jk = 1,jpk
+ IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN
+ tmpnum = tmpnum + (wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk)*tsn(ji,jj,jk,jp_tem))
+ tmpdenom = tmpdenom + (wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk))
+ endif
+ END DO
+ t_zmean(ji,jj) = tmpnum/tmpdenom
+
+ !Depth mean salinity
+ tmpdenom = 0.
+ tmpnum = 0.
+ DO jk = 1,jpk
+ IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN
+ tmpnum = tmpnum + (wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk)*tsn(ji,jj,jk,jp_sal))
+ tmpdenom = tmpdenom + (wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk))
+ endif
+ END DO
+ s_zmean(ji,jj) = tmpnum/tmpdenom
+
+ !save into a 3d grid
+ DO jk = 1,jpk
+ IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN
+ t_zmean_mat(ji,jj,jk) = t_zmean(ji,jj)
+ s_zmean_mat(ji,jj,jk) = s_zmean(ji,jj)
+ else
+ t_zmean_mat(ji,jj,jk) = tsn(ji,jj,jk,jp_tem)
+ s_zmean_mat(ji,jj,jk) = tsn(ji,jj,jk,jp_sal)
+ endif
+ END DO
+ else
+ t_zmean(ji,jj) = tsn(ji,jj,1,jp_tem)
+ s_zmean(ji,jj) = tsn(ji,jj,1,jp_sal)
+ DO jk = 1,jpk
+ t_zmean_mat(ji,jj,jk) = tsn(ji,jj,jk,jp_tem)
+ s_zmean_mat(ji,jj,jk) = tsn(ji,jj,jk,jp_sal)
+ END DO
+ endif
+
+ END DO
+ END DO
+
+ !Calculate the density from the depth varying, and depth average temperature and salinity
+ !-----------------------------
+ !-----------------------------
+
+ ts_pea_mat(:,:,:,:) = tsn(:,:,:,:)
+
+ ts_pea_mat_TS_mean(:,:,:,1) = t_zmean_mat(:,:,:)
+ ts_pea_mat_TS_mean(:,:,:,2) = s_zmean_mat(:,:,:)
+
+ ts_pea_mat_S_mean(:,:,:,1) = t_zmean_mat(:,:,:)
+ ts_pea_mat_S_mean(:,:,:,2) = tsn(:,:,:,jp_sal)
+
+ CALL eos ( ts_pea_mat, tmp_pea_rho, gdept_n(:,:,:) )
+ CALL eos ( ts_pea_mat_TS_mean, tmp_pea_TS_mean_rho, gdept_n(:,:,:) )
+ CALL eos ( ts_pea_mat_S_mean, tmp_pea_S_mean_rho, gdept_n(:,:,:) )
+ tmp_pea_rho = (tmp_pea_rho * rau0) + rau0
+ tmp_pea_TS_mean_rho = (tmp_pea_TS_mean_rho * rau0) + rau0
+ tmp_pea_S_mean_rho = (tmp_pea_S_mean_rho * rau0) + rau0
+
+
+ ! to test the density calculation
+ !CALL iom_put( "tmp_pea_rho" , tmp_pea_rho ) ! pea
+ !CALL iom_put( "tmp_pea_TS_mean_rho" , tmp_pea_TS_mean_rho ) ! pea
+
+
+ ! Caluclation of the PEA.
+ DO jj = 1,jpj
+ DO ji = 1,jpi
+ pea(ji,jj) = 0.
+ peat(ji,jj) = 0.
+ peas(ji,jj) = 0.
+ maxz = 0.
+ int_y_pea = 0.
+ int_y_pea_t = 0.
+ IF ( tmask(ji,jj,1) == 1.0_wp ) THEN ! for sea points
+
+ ! the depth integrated calculation is summed up over the depths, and then divided by the depth
+ DO jk = 1,jpk
+ !for each level...
+
+ IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN ! if above the sea bed...
+ int_y_pea = -((tmp_pea_TS_mean_rho(ji,jj,jk)) - (tmp_pea_rho(ji,jj,jk)))*9.81*gdept_n(ji,jj,jk)*wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk)
+ int_y_pea_t = -((tmp_pea_S_mean_rho(ji,jj,jk)) - (tmp_pea_rho(ji,jj,jk)))*9.81*gdept_n(ji,jj,jk)*wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk)
+ else
+ int_y_pea = 0.
+ int_y_pea_t = 0.
+ endif
+
+ ! check that the sum is not NaN.
+ if ( int_y_pea .ne. int_y_pea ) int_y_pea = 0.
+ if ( int_y_pea_t .ne. int_y_pea_t ) int_y_pea_t = 0.
+ !if ( (int_y_pea*int_y_pea ) .gt. 1.0e6 ) int_y_pea = 0.
+ !if ( (int_y_pea_t*int_y_pea_t) .gt. 1.0e6 ) int_y_pea_t = 0.
+
+ pea(ji,jj) = pea(ji,jj) + int_y_pea
+ peat(ji,jj) = peat(ji,jj) + int_y_pea_t
+ maxz = maxz + (e3t_n(ji,jj,jk)*wgt_co_mat(ji,jj,jk))
+ enddo
+
+
+ !divide by the depth
+ pea(ji,jj) = pea(ji,jj)/maxz
+ peat(ji,jj) = peat(ji,jj)/maxz
+ peas(ji,jj) = pea(ji,jj) - peat(ji,jj)
+
+
+ else
+ pea(ji,jj) = 0.
+ peat(ji,jj) = 0.
+ peas(ji,jj) = 0.
+ endif
+ enddo
+ enddo
+!
+ CALL iom_put( "pea" , pea ) ! pea
+ CALL iom_put( "peat" , peat ) ! pea
+ CALL iom_put( "peas" , peas ) ! pea
+
+
+ DEALLOCATE(ts_pea_mat,ts_pea_mat_TS_mean,ts_pea_mat_S_mean,tmp_pea_rho,tmp_pea_TS_mean_rho,tmp_pea_S_mean_rho)
+
+ END SUBROUTINE dia_pea
+
+END MODULE diapea
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diaregmean.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diaregmean.F90 (revision 10958)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diaregmean.F90 (revision 10958)
@@ -0,0 +1,1103 @@
+MODULE diaregmean
+ !!======================================================================
+ !! *** MODULE diaharm ***
+ !! Timeseries of Regional Means
+ !!======================================================================
+ !! History : 3.6 ! 11/2016 (J Tinker) Original code
+ !!----------------------------------------------------------------------
+ USE oce ! ocean dynamics and tracers variables
+ USE dom_oce ! ocean space and time domain
+ USE in_out_manager ! I/O units
+ USE iom ! I/0 library
+ USE wrk_nemo ! working arrays
+ USE diapea ! PEA
+ USE zdfmxl ! MLD
+ USE sbc_oce
+#if defined key_diaar5
+ USE diaar5
+#endif
+
+
+#if defined key_fabm
+ USE trc
+ USE par_fabm
+#endif
+
+ IMPLICIT NONE
+ PRIVATE
+
+ LOGICAL , PUBLIC :: ln_diaregmean ! region mean calculation
+ PUBLIC dia_regmean_init ! routine called by nemogcm.F90
+ PUBLIC dia_regmean ! routine called by diawri.F90
+ PUBLIC dia_calctmb_region_mean ! routine called by diatmb.F90
+
+
+
+ LOGICAL :: ln_diaregmean_ascii ! region mean calculation ascii output
+ LOGICAL :: ln_diaregmean_bin ! region mean calculation binary output
+ LOGICAL :: ln_diaregmean_nc ! region mean calculation netcdf output
+ LOGICAL :: ln_diaregmean_diaar5 ! region mean calculation including AR5 SLR terms
+ LOGICAL :: ln_diaregmean_diasbc ! region mean calculation including Surface BC
+ LOGICAL :: ln_diaregmean_karamld ! region mean calculation including kara mld terms
+ LOGICAL :: ln_diaregmean_pea ! region mean calculation including pea terms
+
+
+ LOGICAL :: ln_diaregmean_bgc ! region mean calculation including BGC terms
+
+
+
+ REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: tmp_region_mask_real ! tempory region_mask of reals
+ INTEGER, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: region_mask ! region_mask matrix
+ INTEGER :: nmasks ! Number of mask files in region_mask.nc file -
+ INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: nreg_mat ! Number of regions in each mask
+
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: tmp_field_mat !: temporary region_mask
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: tmp_field_AR5_mat !: temporary region_mask
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: tmp_field_SBC_mat !: temporary region_mask
+ INTEGER :: tmp_field_cnt ! tmp_field_cnt integer
+ !!----------------------------------------------------------------------
+ !! NEMO/OPA 3.6 , NEMO Consortium (2014)
+ !! $Id$
+ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+ !!----------------------------------------------------------------------
+CONTAINS
+
+ SUBROUTINE dia_regmean_init
+ !!---------------------------------------------------------------------------
+ !! *** ROUTINE dia_regmean_init ***
+ !!
+ !! ** Purpose: Initialization of region mask namelist
+ !!
+ !! ** Method : Read namelist
+ !! History
+ !! 3.6 ! 11-16 (J Tinker) Routine to initialize dia_regmean
+ !!---------------------------------------------------------------------------
+ !!
+ INTEGER :: ios ! Local integer output status for namelist read
+ INTEGER :: inum ! temporary logical unit ! copied from DOM/domzgr.F90
+ INTEGER :: ierr ! error integer for IOM_get
+ INTEGER :: idmaskvar ! output of iom_varid
+ INTEGER :: maskno ! counter for number of masks
+ INTEGER :: jj,ji ! i and j index
+ INTEGER :: tmpint ! temporary integer
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tmpregion !: temporary region_mask
+ INTEGER, DIMENSION(3) :: zdimsz ! number of elements in each of the 3 dimensions (i.e., lon, lat, no of masks, 297, 375, 4) for an array
+ INTEGER :: zndims ! number of dimensions in an array (i.e. 3, )
+
+
+#if defined key_fabm
+ INTEGER :: js,jl,jn, tmp_dummy
+
+ CHARACTER (len=120) :: tmp_name,tmp_long_name, tmp_unit
+
+ INTEGER :: BGC_nlevs,nBGC_output, bgci
+ CHARACTER(len = 10), ALLOCATABLE, DIMENSION(:) :: BGC_stat_name(:),BGC_lev_name(:),BGC_output_var(:)
+#endif
+
+ !
+ NAMELIST/nam_diaregmean/ ln_diaregmean,ln_diaregmean_ascii,ln_diaregmean_bin,ln_diaregmean_nc,&
+ & ln_diaregmean_karamld, ln_diaregmean_pea,ln_diaregmean_diaar5,ln_diaregmean_diasbc,ln_diaregmean_bgc
+
+
+ ! read in Namelist.
+ !!----------------------------------------------------------------------
+ !
+ REWIND ( numnam_ref ) ! Read Namelist nam_diatmb in referdiatmbence namelist : TMB diagnostics
+ READ ( numnam_ref, nam_diaregmean, IOSTAT=ios, ERR= 901 )
+901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaregmean in reference namelist', lwp )
+
+ REWIND( numnam_cfg ) ! Namelist nam_diatmb in configuration namelist TMB diagnostics
+ READ ( numnam_cfg, nam_diaregmean, IOSTAT = ios, ERR = 902 )
+902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaregmean in configuration namelist', lwp )
+ IF(lwm) WRITE ( numond, nam_diaregmean )
+
+ IF(lwp) THEN ! Control print
+ WRITE(numout,*)
+ WRITE(numout,*) 'dia_regmean_init : Output regional mean Diagnostics'
+ WRITE(numout,*) '~~~~~~~~~~~~'
+ WRITE(numout,*) 'Namelist nam_regmean : set regmeanoutputs '
+ WRITE(numout,*) 'Switch for regmean diagnostics (T) or not (F) ln_diaregmean = ', ln_diaregmean
+ WRITE(numout,*) 'Switch for regmean ascii output (T) or not (F) ln_diaregmean_ascii = ', ln_diaregmean_ascii
+ WRITE(numout,*) 'Switch for regmean binary output (T) or not (F) ln_diaregmean_bin = ', ln_diaregmean_bin
+ WRITE(numout,*) 'Switch for regmean netcdf output (T) or not (F) ln_diaregmean_nc = ', ln_diaregmean_nc
+ WRITE(numout,*) 'Switch for regmean kara mld terms (T) or not (F) ln_diaregmean_karamld = ', ln_diaregmean_karamld
+ WRITE(numout,*) 'Switch for regmean PEA terms (T) or not (F) ln_diaregmean_pea = ', ln_diaregmean_pea
+ WRITE(numout,*) 'Switch for regmean AR5 SLR terms (T) or not (F) ln_diaregmean_diaar5 = ', ln_diaregmean_diaar5
+ WRITE(numout,*) 'Switch for regmean Surface forcing terms (T) or not (F) ln_diaregmean_diasbc = ', ln_diaregmean_diasbc
+ WRITE(numout,*) 'Switch for regmean BioGeoChemistry terms (T) or not (F) ln_diaregmean_bgc = ', ln_diaregmean_bgc
+ ENDIF
+
+
+ ALLOCATE( tmp_field_mat(jpi,jpj,19), STAT= ierr ) !SS/NB/DT/ZA/VA T/S, SSH, MLD, PEA, PEAT, PEAS
+ IF( ierr /= 0 ) CALL ctl_stop( 'tmp_field_mat: failed to allocate tmp_field_mat array' )
+ tmp_field_mat(:,:,:) = 0.
+ tmp_field_cnt = 0
+
+ IF(ln_diaregmean_diaar5) THEN
+ ALLOCATE( tmp_field_AR5_mat(jpi,jpj,4), STAT= ierr ) !SLR terms
+ IF( ierr /= 0 ) CALL ctl_stop( 'tmp_field_AR5_mat: failed to allocate tmp_field_AR5_mat array' )
+ tmp_field_AR5_mat(:,:,:) = 0.
+ ENDIF
+
+ IF(ln_diaregmean_diasbc) THEN
+ ALLOCATE( tmp_field_SBC_mat(jpi,jpj,7), STAT= ierr ) !SBC terms
+ IF( ierr /= 0 ) CALL ctl_stop( 'tmp_field_SBC_mat: failed to allocate tmp_field_SBC_mat array' )
+ tmp_field_SBC_mat(:,:,:) = 0.
+ ENDIF
+
+
+#if defined key_fabm
+ ! as there are so many BGC variables, write out the necessary iodef.xml and field_def.xml entries into ocean.output
+
+ IF(ln_diaregmean_bgc) THEN
+ IF(lwp) THEN ! Control print
+
+ BGC_nlevs = 5
+ ALLOCATE( BGC_stat_name(6),BGC_lev_name(BGC_nlevs))
+ nBGC_output = 16
+ ALLOCATE( BGC_output_var(nBGC_output))
+
+ BGC_output_var(1) = 'N1_p'
+ BGC_output_var(2) = 'N3_n'
+ BGC_output_var(3) = 'N4_n'
+ BGC_output_var(4) = 'N5_s'
+ BGC_output_var(5) = 'O2_o'
+ BGC_output_var(6) = 'P1_Chl'
+ BGC_output_var(7) = 'P2_Chl'
+ BGC_output_var(8) = 'P3_Chl'
+ BGC_output_var(9) = 'P4_Chl'
+ BGC_output_var(10) = 'P1_c'
+ BGC_output_var(11) = 'P2_c'
+ BGC_output_var(12) = 'P3_c'
+ BGC_output_var(13) = 'P4_c'
+ BGC_output_var(14) = 'Z4_c'
+ BGC_output_var(15) = 'Z5_c'
+ BGC_output_var(16) = 'Z6_c'
+
+ BGC_stat_name(1) = '_ave'
+ BGC_stat_name(2) = '_tot'
+ BGC_stat_name(3) = '_var'
+ BGC_stat_name(4) = '_cnt'
+ BGC_stat_name(5) = '_reg_id'
+ BGC_stat_name(6) = '_mask_id'
+ BGC_lev_name(1) = 'top'
+ BGC_lev_name(2) = 'bot'
+ BGC_lev_name(3) = 'dif'
+ BGC_lev_name(4) = 'zav'
+ BGC_lev_name(5) = 'vol'
+
+
+ WRITE(numout,*) ''
+ WRITE(numout,*) 'diaregmean BGC field_def.xml entries'
+ WRITE(numout,*) ''
+
+
+ DO jn=1,jp_fabm ! State loop
+ DO js=1,6
+ DO jl=1,BGC_nlevs
+
+ tmp_name=TRIM( TRIM("reg_")//TRIM(BGC_lev_name(jl))//TRIM("_")//TRIM(ctrcnm(jn))// TRIM(BGC_stat_name(js)) )
+
+ tmp_long_name = TRIM(ctrcln(jn))
+ tmp_unit = TRIM(ctrcun(jn))
+
+ ! Where using volume integrated values, change units...
+
+ IF ((jl .EQ. 5) .AND. (js .EQ. 2)) then
+ SELECT CASE (trim(tmp_unit))
+ CASE ('mg C/m^3') ; tmp_unit = 'Mg C (T C)' !'mg C/m^3'
+ CASE ('mg/m^3') ; tmp_unit = 'Mg (T)' !'mg/m^3'
+ CASE ('mmol C/m^3') ; tmp_unit = 'Mmol C' !'mmol C/m^3'
+ CASE ('mmol N/m^3') ; tmp_unit = 'Mmol N' !'mmol N/m^3'
+ CASE ('mmol O_2/m^3') ; tmp_unit = 'Mmol O' !'mmol O_2/m^3'
+ CASE ('mmol P/m^3') ; tmp_unit = 'Mmol P' !'mmol P/m^3'
+ CASE ('mmol Si/m^3') ; tmp_unit = 'Mmol S' !'mmol Si/m^3'
+ CASE ('umol/kg') ; tmp_unit = 'Mmol' !'umol/kg' = mmol/m^3
+ ! CASE ('1/m') ; cycle
+ CASE DEFAULT
+ tmp_unit = TRIM(TRIM(tmp_unit)//TRIM('x 1e9 m^3'))
+ END SELECT
+ ENDIF
+
+ WRITE(numout,*) TRIM(TRIM(''))
+
+ END DO
+ END DO
+ END DO
+
+ WRITE(numout,*) ''
+ WRITE(numout,*) 'diaregmean BGC iodef.xml entries'
+ WRITE(numout,*) ''
+ DO js=1,6
+
+ DO jn=1,jp_fabm ! State loop
+
+ DO bgci=1,nBGC_output!
+ if (trim(ctrcnm(jn)) == TRIM(BGC_output_var(bgci))) CYCLE
+ ENDDO
+ DO jl=1,BGC_nlevs
+ ! only print out area averages for ss, nb, diff, and depth averaged, and total values for volume integrated
+ IF ((jl .EQ. 5) .AND. (js .NE. 2)) CYCLE ! cycle if vol, and not tot.
+ IF ((jl .NE. 5) .AND. (js .NE. 1)) CYCLE ! cycle if other levels, and not ave.
+
+ tmp_name=TRIM(TRIM("reg_")//TRIM(BGC_lev_name(jl))//TRIM("_")//TRIM(ctrcnm(jn))// TRIM(BGC_stat_name(js)))
+ tmp_long_name = TRIM(ctrcln(jn))
+
+ WRITE(numout,*) TRIM(TRIM(''))
+
+ END DO !level
+ END DO ! State loop
+ END DO !statistic
+ WRITE(numout,*) ''
+ DEALLOCATE( BGC_stat_name,BGC_lev_name)
+
+ ENDIF ! Control print
+
+ ENDIF !ln_diaregmean_bgc
+
+#endif
+
+
+ IF (ln_diaregmean) THEN
+
+ ! Open region mask for region means, and retrieve the size of the mask (number of levels)
+ CALL iom_open ( 'region_mask.nc', inum )
+ idmaskvar = iom_varid( inum, 'mask', kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE.)
+ nmasks = zdimsz(3)
+
+ ! read in the region mask (which contains floating point numbers) into a temporary array of reals.
+ ALLOCATE( tmp_region_mask_real(jpi,jpj,nmasks), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_regmean_init: failed to allocate tmp_region_mask_real array' )
+
+ ! Use jpdom_unknown to read in a n-layer mask.
+ tmp_region_mask_real(:,:,:) = 0
+ CALL iom_get( inum, jpdom_unknown, 'mask', tmp_region_mask_real(1:nlci,1:nlcj,1:nmasks), &
+ & kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nmasks /) )
+
+ CALL iom_close( inum )
+
+ !Convert the region mask of reals into one of integers.
+
+ ALLOCATE( region_mask(jpi,jpj,nmasks), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_regmean_init: failed to allocate region_mask array' )
+ region_mask(:,:,:) = 0
+ region_mask = int(tmp_region_mask_real(:,:,:))
+ DEALLOCATE( tmp_region_mask_real)
+
+
+ ALLOCATE( nreg_mat(nmasks), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_regmean_init: failed to allocate nreg_mat array' )
+
+ ! work out the number of regions in each mask, asssuming land is 0, and the regions are consectively numbered,
+ ! without missing any number, so the number of regions is the maximum number + 1 (for land). mpp_max across the
+ ! processors to get the global maxima
+ DO maskno = 1,nmasks
+ tmpint = maxval(region_mask(:,:,maskno))
+ CALL mpp_max( tmpint )
+ nreg_mat(maskno) = tmpint + 1
+ END DO
+
+ IF(lwp) THEN
+ ! if writing out as binary and text, open the files.
+ IF ( ln_diaregmean_bin ) THEN
+ ! Open binary for region means
+ CALL ctl_opn( numdct_reg_bin ,'region_mean_timeseries.dat' , 'NEW', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. )
+ ENDIF
+
+ IF ( ln_diaregmean_ascii ) THEN
+ ! Open text files for region means
+ CALL ctl_opn( numdct_reg_txt ,'region_mean_timeseries.txt' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. )
+ ENDIF
+ ENDIF
+ ENDIF
+
+ END SUBROUTINE dia_regmean_init
+
+ SUBROUTINE dia_calctmb_region_mean( pinfield,pouttmb )
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE dia_calctmb_region_mean ***
+ !!
+ !! ** Purpose : Find the Top, Bottom and Top minus Bottom fields of water Column
+ !! : and depth average, and volume and mass intergated values.
+
+ !!
+ !! ** Method :
+ !! use mbathy to find surface, mid and bottom of model levels
+ !!
+ !! History :
+ !! 3.6 ! 08-14 (E. O'Dea) Routine based on dia_wri_foam
+ !!----------------------------------------------------------------------
+ !! * Modules used
+
+ ! Routine to map 3d field to top, middle, bottom
+ IMPLICIT NONE
+
+
+ ! Routine arguments
+ REAL(wp), DIMENSION(jpi, jpj, jpk), INTENT(IN ) :: pinfield ! Input 3d field and mask
+ REAL(wp), DIMENSION(jpi, jpj, 6 ), INTENT( OUT) :: pouttmb ! Output top, bottom and surface minus bed, zav, vol int, mass int
+
+ ! Local variables
+ INTEGER :: ji,jj,jk ! Dummy loop indices
+
+ ! Local Real
+ REAL(wp) :: zmdi ! set masked values
+ ! for depth int
+ REAL(wp) :: tmpnumer,tmpnumer_mass,tmpdenom ,z_av_val,vol_int_val
+
+ zmdi=1.e+20 !missing data indicator for masking
+
+ !zmdi=0 !missing data indicator for masking
+
+ ! Calculate top
+ pouttmb(:,:,1) = pinfield(:,:,1)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+
+ ! Calculate middle
+ !DO jj = 1,jpj
+ ! DO ji = 1,jpi
+ ! jk = max(1,mbathy(ji,jj)/2)
+ ! pouttmb(ji,jj,2) = pinfield(ji,jj,jk)*tmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+ ! END DO
+ !END DO
+
+ ! Calculate bottom, and top minus bottom
+ DO jj = 1,jpj
+ DO ji = 1,jpi
+ IF ( tmask(ji,jj,1) .EQ. 1) THEN ! if land
+
+ jk = max(1,mbathy(ji,jj) - 1)
+ pouttmb(ji,jj,2) = pinfield(ji,jj,jk)*tmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+
+ pouttmb(ji,jj,3) = (pouttmb(ji,jj,1) - pouttmb(ji,jj,2))*tmask(ji,jj,1) + zmdi*(1.0-tmask(ji,jj,1))
+
+ !Depth and volume integral:
+ !---------------------------
+ !Vol int = Concentration * vol of grid box, summed over depth.
+ !Mass int = Concentration * vol of grid box * density of water, summed over depth.
+ !Depth Average = Vol int divided by * (vol of grid box summed over depth).
+
+ tmpnumer = 0.
+ tmpnumer_mass = 0.
+ tmpdenom = 0.
+ DO jk = 1,jpk
+ tmpnumer = tmpnumer + pinfield(ji,jj,jk)*tmask(ji,jj,jk)*e1t(ji,jj)*e2t(ji,jj)*e3t_n(ji,jj,jk)
+ tmpnumer_mass = tmpnumer_mass + pinfield(ji,jj,jk)*tmask(ji,jj,jk)*e1t(ji,jj)*e2t(ji,jj)*e3t_n(ji,jj,jk)*rhop(ji,jj,jk)
+ tmpdenom = tmpdenom + tmask(ji,jj,jk)*e1t(ji,jj)*e2t(ji,jj)*e3t_n(ji,jj,jk)
+ END DO
+ !z_av_val = tmpnumer/tmpdenom
+ !vol_int_val = tmpnumer
+ !mass_int_val = tmpnumer*density
+
+ pouttmb(ji,jj,4) = tmpnumer/tmpdenom ! depth averaged
+ pouttmb(ji,jj,5) = tmpnumer ! Vol integrated
+ pouttmb(ji,jj,6) = tmpnumer_mass ! Mass integrated (for heat and salt calcs)
+ ELSE
+ pouttmb(ji,jj,1) = zmdi
+ pouttmb(ji,jj,2) = zmdi
+ pouttmb(ji,jj,3) = zmdi
+ pouttmb(ji,jj,4) = zmdi
+ pouttmb(ji,jj,5) = zmdi
+ pouttmb(ji,jj,6) = zmdi
+ ENDIF
+ END DO
+ END DO
+
+ END SUBROUTINE dia_calctmb_region_mean
+
+
+ SUBROUTINE dia_regmean( kt )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE dia_regmean ***
+ !! ** Purpose : Produce regional mean diagnostics
+ !!
+ !! ** Method : calls dia_wri_region_mean to calculate and write the regional means for a number of variables,
+ !! (calling dia_calctmb_region_mean where necessary).
+ !!
+ !! Closes all text and binary files on last time step
+ !!
+ !!
+ !!
+ !!
+ !! History :
+ !! 3.6 ! 11-16 (J. Tinker)
+ !!
+ !!--------------------------------------------------------------------
+ REAL(wp), POINTER, DIMENSION(:,:,:) :: tmp1mat ! temporary array of 1's
+ REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmbT ! temporary T workspace
+ REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmbS ! temporary S workspace
+ REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmb1 ! temporary density workspace
+ REAL(wp) :: zmdi ! set masked values
+ INTEGER, INTENT( in ) :: kt ! ocean time-step index
+
+ REAL(wp) :: zdt ! temporary reals
+ INTEGER :: i_steps, ierr ! no of timesteps per hour, allocation error index
+ INTEGER :: maskno,jj,ji,jk,jm,nreg ! indices of mask, i and j, and number of regions
+
+#if defined key_fabm
+ INTEGER :: jn ,tmp_dummy ! set masked values
+ REAL(wp) :: tmp_val ! tmp value, to allow min and max value clamping (not implemented)
+ INTEGER :: jl
+ CHARACTER (len=60) :: tmp_name_bgc_top,tmp_name_bgc_bot,tmp_name_bgc_dif, tmp_name_bgc_zav, tmp_name_bgc_vol
+ CHARACTER (len=60) :: tmp_output_filename
+ REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmbBGC ! temporary BGC workspace
+
+ LOGICAL :: verbose
+ verbose = .FALSE.
+ tmp_val = 0
+#endif
+ zmdi=1.e+20 !missing data indicator for maskin
+
+ IF (ln_diaregmean) THEN
+ ! If regional mean calculations required by namelist
+ ! -----------------
+ ! identify hourly time steps (not used)
+ zdt = rdt
+ IF( nacc == 1 ) zdt = rdtmin
+
+ IF( MOD( 3600,INT(zdt) ) == 0 ) THEN
+ i_steps = 3600/INT(zdt)
+ ELSE
+ CALL ctl_stop('STOP', 'dia_regmean: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible')
+ ENDIF
+
+ ! Every time step, add physical, SBC, PEA, MLD terms to create hourly sums.
+ ! Every hour, then hourly sums are divided by the number of timesteps in the hour to make hourly means
+ ! These hourly mean values are then used to caluclate the regional means, and output with IOM.
+#if defined key_fabm
+ ! BGC values are not averaged up over the hour, but are output as hourly instantaneous values.
+#endif
+
+
+ !Extract 2d fields from 3d T and S with dia_calctmb_region_mean
+ CALL wrk_alloc( jpi , jpj, 6 , zwtmbT )
+ CALL wrk_alloc( jpi , jpj, 6 , zwtmbS )
+ CALL wrk_alloc( jpi , jpj, 6 , zwtmb1 )
+
+ CALL dia_calctmb_region_mean( tsn(:,:,:,jp_tem),zwtmbT)
+ CALL dia_calctmb_region_mean( tsn(:,:,:,jp_sal),zwtmbS)
+
+ ! To calc regional mean time series of int vol and mass, run region mean code on array of 1's...
+ ! - then when multplying by volume, gives volume,
+ ! - then when multplying by volume*density, gives mass
+
+ CALL wrk_alloc( jpi , jpj, jpk , tmp1mat )
+ DO jj = 1,jpj
+ DO ji = 1,jpi
+ DO jk = 1,jpk
+ tmp1mat(ji,jj,jk) = 1
+ END DO
+ END DO
+ END DO
+
+ CALL dia_calctmb_region_mean( tmp1mat,zwtmb1)
+ CALL wrk_dealloc( jpi , jpj, jpk , tmp1mat )
+
+ ! Add 2d fields every time step to the hourly total.
+
+ tmp_field_mat(:,:,1) = tmp_field_mat(:,:,1) + (zwtmbT(:,:,1)*tmask(:,:,1)) !sst
+ tmp_field_mat(:,:,2) = tmp_field_mat(:,:,2) + (zwtmbT(:,:,2)*tmask(:,:,1)) !nbt
+ tmp_field_mat(:,:,3) = tmp_field_mat(:,:,3) + (zwtmbT(:,:,3)*tmask(:,:,1)) !dft
+
+ tmp_field_mat(:,:,4) = tmp_field_mat(:,:,4) + (zwtmbT(:,:,4)*tmask(:,:,1)) !zat
+ tmp_field_mat(:,:,5) = tmp_field_mat(:,:,5) + (zwtmbT(:,:,5)*tmask(:,:,1)) !vat
+ tmp_field_mat(:,:,6) = tmp_field_mat(:,:,6) + ((zwtmbT(:,:,6)*tmask(:,:,1)*4.2e3))! heat
+
+ tmp_field_mat(:,:,7) = tmp_field_mat(:,:,7) + (zwtmbS(:,:,1)*tmask(:,:,1)) !sss
+ tmp_field_mat(:,:,8) = tmp_field_mat(:,:,8) + (zwtmbS(:,:,2)*tmask(:,:,1)) !nbs
+ tmp_field_mat(:,:,9) = tmp_field_mat(:,:,9) + (zwtmbS(:,:,3)*tmask(:,:,1)) !dfs
+
+ tmp_field_mat(:,:,10) = tmp_field_mat(:,:,10) + (zwtmbS(:,:,4)*tmask(:,:,1)) !zas
+ tmp_field_mat(:,:,11) = tmp_field_mat(:,:,11) + (zwtmbS(:,:,5)*tmask(:,:,1)) !vas
+ tmp_field_mat(:,:,12) = tmp_field_mat(:,:,12) + (zwtmbS(:,:,6)*tmask(:,:,1)) !salt
+
+ tmp_field_mat(:,:,13) = tmp_field_mat(:,:,13) + (zwtmb1(:,:,5)*tmask(:,:,1))!vol
+ tmp_field_mat(:,:,14) = tmp_field_mat(:,:,14) + (zwtmb1(:,:,6)*tmask(:,:,1))!mass
+
+ tmp_field_mat(:,:,15) = tmp_field_mat(:,:,15) + (sshn(:,:)*tmask(:,:,1)) !ssh
+
+ CALL wrk_dealloc( jpi , jpj, 6 , zwtmbT )
+ CALL wrk_dealloc( jpi , jpj, 6 , zwtmbS )
+ CALL wrk_dealloc( jpi , jpj, 6 , zwtmb1 )
+
+ IF( ln_diaregmean_karamld ) THEN
+ tmp_field_mat(:,:,16) = tmp_field_mat(:,:,16) + (hmld_kara(:,:)*tmask(:,:,1)) !mldkara
+ ENDIF
+
+ IF( ln_diaregmean_pea ) THEN
+ tmp_field_mat(:,:,17) = tmp_field_mat(:,:,17) + (pea(:,:)*tmask(:,:,1)) !pea
+ tmp_field_mat(:,:,18) = tmp_field_mat(:,:,18) + (peat(:,:)*tmask(:,:,1)) !peat
+ tmp_field_mat(:,:,19) = tmp_field_mat(:,:,19) + (peas(:,:)*tmask(:,:,1)) !peas
+ ENDIF
+
+ IF( ln_diaregmean_diaar5 ) THEN
+ tmp_field_AR5_mat(:,:,1) = tmp_field_AR5_mat(:,:,1) + (sshsteric_mat(:,:)*tmask(:,:,1))
+ tmp_field_AR5_mat(:,:,2) = tmp_field_AR5_mat(:,:,2) + (sshthster_mat(:,:)*tmask(:,:,1))
+ tmp_field_AR5_mat(:,:,3) = tmp_field_AR5_mat(:,:,3) + (sshhlster_mat(:,:)*tmask(:,:,1))
+ tmp_field_AR5_mat(:,:,4) = tmp_field_AR5_mat(:,:,4) + (zbotpres_mat(:,:)*tmask(:,:,1))
+ ENDIF
+
+ IF( ln_diaregmean_diasbc ) THEN
+ tmp_field_SBC_mat(:,:,1) = tmp_field_SBC_mat(:,:,1) + ((qsr + qns)*tmask(:,:,1))
+ tmp_field_SBC_mat(:,:,2) = tmp_field_SBC_mat(:,:,2) + (qsr*tmask(:,:,1))
+ tmp_field_SBC_mat(:,:,3) = tmp_field_SBC_mat(:,:,3) + (qns*tmask(:,:,1))
+ tmp_field_SBC_mat(:,:,4) = tmp_field_SBC_mat(:,:,4) + (emp*tmask(:,:,1))
+ tmp_field_SBC_mat(:,:,5) = tmp_field_SBC_mat(:,:,5) + (wndm*tmask(:,:,1))
+ tmp_field_SBC_mat(:,:,6) = tmp_field_SBC_mat(:,:,6) + (pressnow*tmask(:,:,1))
+ tmp_field_SBC_mat(:,:,7) = tmp_field_SBC_mat(:,:,7) + (rnf*tmask(:,:,1))
+ ENDIF
+
+
+
+ tmp_field_cnt = tmp_field_cnt + 1
+
+ ! On the hour, calculate hourly means from the hourly total,and process the regional means.
+
+ IF ( MOD( kt, i_steps ) == 0 .and. kt .ne. nn_it000 ) THEN
+
+
+ CALL dia_wri_region_mean(kt, "sst" , tmp_field_mat(:,:,1)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "nbt" , tmp_field_mat(:,:,2)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "dft" , tmp_field_mat(:,:,3)/real(tmp_field_cnt,wp))
+
+ CALL dia_wri_region_mean(kt, "zat" , tmp_field_mat(:,:,4)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "vat" , tmp_field_mat(:,:,5)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "heat" , tmp_field_mat(:,:,6)/real(tmp_field_cnt,wp)/1e12)
+
+ CALL dia_wri_region_mean(kt, "sss" , tmp_field_mat(:,:,7)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "nbs" , tmp_field_mat(:,:,8)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "dfs" , tmp_field_mat(:,:,9)/real(tmp_field_cnt,wp))
+
+ CALL dia_wri_region_mean(kt, "zas" , tmp_field_mat(:,:,10)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "vas" , tmp_field_mat(:,:,11)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "salt" , tmp_field_mat(:,:,12)/real(tmp_field_cnt,wp)/1e12)
+
+ CALL dia_wri_region_mean(kt, "vol" , tmp_field_mat(:,:,13)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "mass" , tmp_field_mat(:,:,14)/real(tmp_field_cnt,wp))
+
+ CALL dia_wri_region_mean(kt, "ssh" , tmp_field_mat(:,:,15)/real(tmp_field_cnt,wp))
+
+
+ IF( ln_diaregmean_karamld ) THEN
+ CALL dia_wri_region_mean(kt, "mldkara" , tmp_field_mat(:,:,16)/real(tmp_field_cnt,wp)) ! tm
+ ENDIF
+
+ IF( ln_diaregmean_pea ) THEN
+ CALL dia_wri_region_mean(kt, "pea" , tmp_field_mat(:,:,17)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "peat" , tmp_field_mat(:,:,18)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "peas" , tmp_field_mat(:,:,19)/real(tmp_field_cnt,wp)) ! tmb
+ ENDIF
+
+ tmp_field_mat(:,:,:) = 0.
+
+ IF( ln_diaregmean_diaar5 ) THEN
+
+ CALL dia_wri_region_mean(kt, "ssh_steric" , tmp_field_AR5_mat(:,:,1)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "ssh_thermosteric", tmp_field_AR5_mat(:,:,2)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "ssh_halosteric" , tmp_field_AR5_mat(:,:,3)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "bot_pres" , tmp_field_AR5_mat(:,:,4)/real(tmp_field_cnt,wp))
+ tmp_field_AR5_mat(:,:,:) = 0.
+ ENDIF
+
+ IF( ln_diaregmean_diasbc ) THEN
+
+ CALL dia_wri_region_mean(kt, "qt" , tmp_field_SBC_mat(:,:,1)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "qsr" , tmp_field_SBC_mat(:,:,2)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "qns" , tmp_field_SBC_mat(:,:,3)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "emp" , tmp_field_SBC_mat(:,:,4)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "wspd" , tmp_field_SBC_mat(:,:,5)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "mslp" , tmp_field_SBC_mat(:,:,6)/real(tmp_field_cnt,wp))
+ CALL dia_wri_region_mean(kt, "rnf" , tmp_field_SBC_mat(:,:,7)/real(tmp_field_cnt,wp))
+ tmp_field_SBC_mat(:,:,:) = 0.
+ ENDIF
+
+#if defined key_fabm
+ !ADD Biogeochemistry
+
+ IF( ln_diaregmean_bgc ) THEN !ln_diaregmean_bgc
+
+ ! Loop through 3d BGC tracers
+ DO jn=1,jp_fabm ! State loop
+
+ ! get variable name for different levels
+ tmp_name_bgc_top=TRIM(TRIM("top_")//TRIM(ctrcnm(jn)))
+ tmp_name_bgc_bot=TRIM(TRIM("bot_")//TRIM(ctrcnm(jn)))
+ tmp_name_bgc_dif=TRIM(TRIM("dif_")//TRIM(ctrcnm(jn)))
+ tmp_name_bgc_zav=TRIM(TRIM("zav_")//TRIM(ctrcnm(jn)))
+ tmp_name_bgc_vol=TRIM(TRIM("vol_")//TRIM(ctrcnm(jn)))
+
+ ! print out names if verbose
+ IF(verbose .AND. lwp) THEN
+ WRITE(numout,*)
+ WRITE(numout,*) 'dia_regmean tmp_name_bgc_top : ',TRIM(tmp_name_bgc_top)
+ WRITE(numout,*) 'dia_regmean tmp_name_bgc_bot : ',TRIM(tmp_name_bgc_bot)
+ WRITE(numout,*) 'dia_regmean tmp_name_bgc_dif : ',TRIM(tmp_name_bgc_dif)
+ WRITE(numout,*) 'dia_regmean tmp_name_bgc_zav : ',TRIM(tmp_name_bgc_zav)
+ WRITE(numout,*) 'dia_regmean tmp_name_bgc_vol : ',TRIM(tmp_name_bgc_vol)
+ CALL FLUSH(numout)
+
+ ENDIF
+
+ !Allocate working array, and get surface, bed etc fields.
+ CALL wrk_alloc( jpi , jpj, 6 , zwtmbBGC )
+ CALL dia_calctmb_region_mean( trn(:,:,:,jn),zwtmbBGC )
+
+
+ !Print out 2d fields to ascii text files to check values if verbose. (24MB per time step, per BGC variable)
+ IF (verbose) THEN
+
+ WRITE (tmp_output_filename, "(A4,I3.3,A1,I6.6,A1,I3.3,A4)") "bgc_",jn,"_",kt,"_",narea,".txt"
+ WRITE (*,*) tmp_output_filename
+ OPEN(UNIT=74,FILE=TRIM(tmp_output_filename))
+
+ DO ji = 1,jpi
+ DO jj = 1,jpj
+ WRITE(74,FMT="(I4,I4,F3,F25.5,F25.5,F25.5,F25.5,F25.5)") nimpp+ji, njmpp+jj,tmask(ji,jj,1),&
+ & zwtmbBGC(ji,jj,1),zwtmbBGC(ji,jj,2),zwtmbBGC(ji,jj,3),zwtmbBGC(ji,jj,4),zwtmbBGC(ji,jj,5)/1e9
+ END DO
+ END DO
+ CLOSE(74)
+ ENDIF
+
+ ! Do region means
+ CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_top) , zwtmbBGC(:,:,1))
+ CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_bot) , zwtmbBGC(:,:,2))
+ CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_dif) , zwtmbBGC(:,:,3))
+ CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_zav) , zwtmbBGC(:,:,4))
+ CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_vol) , zwtmbBGC(:,:,5)/1e9)
+
+
+ !Deallocate working array
+ CALL wrk_dealloc( jpi , jpj, 6 , zwtmbBGC )
+ ENDDO ! State loop
+ ENDIF !ln_diaregmean_bgc
+
+#endif
+
+ tmp_field_cnt = 0
+
+ ENDIF ! ( MOD( kt, i_steps ) == 0 .and. kt .ne. nn_it000 )
+
+
+ ! If on the last time step, close binary and ascii files.
+ IF( kt == nitend ) THEN
+ IF(lwp) THEN
+ IF ( ln_diaregmean_bin ) THEN
+ !Closing binary files for regional mean time series.
+ CLOSE(numdct_reg_bin)
+ ENDIF
+ IF ( ln_diaregmean_ascii ) THEN
+ !Closing text files for regional mean time series.
+ CLOSE(numdct_reg_txt)
+ ENDIF
+
+ DEALLOCATE( region_mask, nreg_mat, tmp_field_mat)
+ IF( ln_diaregmean_diaar5 ) DEALLOCATE( tmp_field_AR5_mat)
+ IF( ln_diaregmean_diasbc ) DEALLOCATE( tmp_field_SBC_mat)
+ ENDIF
+ ENDIF
+
+
+ ELSE
+ CALL ctl_warn('dia_regmean: regmean diagnostic is set to false you should not have seen this')
+ ENDIF
+
+ END SUBROUTINE dia_regmean
+
+
+ SUBROUTINE dia_wri_region_mean(kt, tmp_name, infield )
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE dia_tmb ***
+ !!
+ !! ** Purpose : Calculate and write region mean time series for 2d arrays
+ !!
+ !! ** Method :
+ !! use
+ !!
+ !! History :
+ !! ?? ! 15/10/2015 (JTinker) Routine taken from old dia_wri_foam
+ !!----------------------------------------------------------------------
+ !! * Modules used
+ !use lib_mpp
+ !use lib_fortr
+ IMPLICIT NONE
+
+ INTEGER, INTENT(in) :: kt
+ CHARACTER (len=*) , INTENT(IN ) :: tmp_name
+ REAL(wp), DIMENSION(jpi, jpj), INTENT(IN ) :: infield ! Input 3d field and mask
+
+ ! Local variables
+ INTEGER, DIMENSION(jpi, jpj) :: internal_region_mask ! Input 3d field and mask
+ REAL(wp), DIMENSION(jpi, jpj) :: internal_infield ! Internal data field
+ REAL(wp), ALLOCATABLE, DIMENSION(:) :: zrmet_ave,zrmet_tot,zrmet_var,zrmet_cnt,zrmet_mask_id,zrmet_reg_id ,zrmet_min,zrmet_max
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zrmet_out
+ REAL(wp), ALLOCATABLE, DIMENSION(:) :: ave_mat,tot_mat,num_mat,var_mat,ssq_mat,cnt_mat,reg_id_mat,mask_id_mat ,min_mat,max_mat !: region_mask
+
+ REAL(wp) :: zmdi, zrmet_val ! set masked values
+ INTEGER :: maskno,nreg ! ocean time-step indexocean time step
+ INTEGER :: ji,jj,jk,ind,jm ! Dummy loop indices
+ INTEGER :: reg_ind_cnt ! Dummy loop indices
+
+ INTEGER :: ierr
+ REAL(wp) :: tmpreal
+ CHARACTER(LEN=180) :: FormatString,nreg_string,tmp_name_iom
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: dummy_zrmet
+ LOGICAL :: verbose
+ verbose = .False.
+
+
+ zmdi=1.e+20 !missing data indicator for maskin
+
+ !Allocate output arrays for iomput, set to zmdi, and set a region counter = 1
+ ALLOCATE( zrmet_ave(n_regions_output), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_ave array' )
+ ALLOCATE( zrmet_tot(n_regions_output), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_tot array' )
+ ALLOCATE( zrmet_var(n_regions_output), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_var array' )
+ ALLOCATE( zrmet_cnt(n_regions_output), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_cnt array' )
+ ALLOCATE( zrmet_mask_id(n_regions_output), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_mask_id array' )
+ ALLOCATE( zrmet_reg_id(n_regions_output), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_reg_id array' )
+
+
+ ALLOCATE( zrmet_min(n_regions_output), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_min array' )
+ ALLOCATE( zrmet_max(n_regions_output), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_max array' )
+
+ ALLOCATE( zrmet_out(jpi,jpj,n_regions_output), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_reg_id array' )
+
+
+
+ IF(lwp .AND. verbose) THEN
+ WRITE(numout,*)
+ WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//';'
+ WRITE(numout,*)
+ ENDIF
+
+ DO ji = 1,jpi
+ DO jj = 1,jpj
+ internal_infield(ji,jj) = infield(ji,jj)
+ END DO
+ END DO
+
+ ! Check for NANS # JT 03/09/2018
+ DO ji = 1,jpi
+ DO jj = 1,jpj
+ IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
+ IF ( internal_infield(ji,jj) .ne. internal_infield(ji,jj) ) THEN
+ WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//' Nan at (kt,i,j): ',kt,ji - (-jpizoom+1-nimpp+1),jj - (-jpjzoom+1-njmpp+1)
+ internal_infield(ji,jj) = 0.
+ ENDIF
+ ELSE
+ IF ( internal_infield(ji,jj) .ne. internal_infield(ji,jj) ) THEN
+ WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//' Masked Nan at (kt,i,j): ',kt,ji - (-jpizoom+1-nimpp+1),jj - (-jpjzoom+1-njmpp+1)
+ internal_infield(ji,jj) = 0.
+ ENDIF
+ ENDIF
+ END DO
+ END DO
+
+
+ zrmet_ave(:) = zmdi
+ zrmet_tot(:) = zmdi
+ zrmet_var(:) = zmdi
+ zrmet_cnt(:) = zmdi
+ zrmet_mask_id(:) = zmdi
+ zrmet_reg_id(:) = zmdi
+
+ zrmet_min(:) = zmdi
+ zrmet_max(:) = zmdi
+ reg_ind_cnt = 1
+
+
+ ! loop though the masks
+ DO maskno = 1,nmasks
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; begin mask loops: ',maskno
+
+
+ ! For each mask, get the number of regions (nreg), and a local copy of the region.
+ nreg = nreg_mat(maskno)
+ internal_region_mask = region_mask(:,:,maskno)
+
+ ! allocate temporary stat arrays, and set to zero
+ ALLOCATE( ave_mat(nreg), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate ave_mat array' )
+ ALLOCATE( tot_mat(nreg), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate tot_mat array' )
+ ALLOCATE( num_mat(nreg), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate num_mat array' )
+ ALLOCATE( var_mat(nreg), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate var_mat array' )
+ ALLOCATE( ssq_mat(nreg), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate ssq_mat array' )
+ ALLOCATE( cnt_mat(nreg), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate cnt_mat array' )
+
+ ALLOCATE( min_mat(nreg), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate min_mat array' )
+ ALLOCATE( max_mat(nreg), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate max_mat array' )
+
+ ALLOCATE( reg_id_mat(nreg), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate reg_id_mat array' )
+ ALLOCATE( mask_id_mat(nreg), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate mask_id_mat array' )
+
+
+
+ ave_mat(:) = 0.
+ tot_mat(:) = 0.
+ num_mat(:) = 0.
+ var_mat(:) = 0.
+ cnt_mat(:) = 0.
+ ssq_mat(:) = 0.
+
+ min_mat(:) = zmdi
+ max_mat(:) = -zmdi
+ reg_id_mat(:) = 0.
+ mask_id_mat(:) = 0.
+
+ ! loop though the array. for each sea grid box where tmask == 1),
+ ! read which region the grid box is in, add the value of the gridbox (and its square)
+ ! to the total for that region, and then increment the counter for that region.
+ !CALL cpu_time(start_reg_mean_loop)
+ !WRITE(numout,*) kt,start_reg_mean_loop
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; begin spatial loops: '
+ DO ji = nldi,nlei
+ DO jj = nldj,nlej
+ IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
+ ind = internal_region_mask(ji,jj)+1
+ tot_mat(ind) = tot_mat(ind) + (internal_infield(ji,jj))
+ ssq_mat(ind) = ssq_mat(ind) + ( internal_infield(ji,jj) * internal_infield(ji,jj))
+ cnt_mat(ind) = cnt_mat(ind) + 1.
+
+ min_mat(ind) = min(min_mat(ind),internal_infield(ji,jj))
+ max_mat(ind) = max(max_mat(ind),internal_infield(ji,jj))
+ ENDIF
+ END DO
+ END DO
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finish spatial loops: '
+ ! sum the totals, the counts, and the squares across the processors
+ CALL mpp_sum( tot_mat,nreg )
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum 1'
+ CALL mpp_sum( ssq_mat,nreg )
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum 2'
+ CALL mpp_sum( cnt_mat,nreg )
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum 2'
+
+ CALL mpp_min( min_mat,nreg )
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_min'
+ CALL mpp_max( max_mat,nreg )
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_max'
+
+
+ !calculate the mean and variance from the total, sum of squares and the count.
+
+ ave_mat = tot_mat(:)/cnt_mat(:)
+ var_mat = ssq_mat(:)/cnt_mat(:) - (ave_mat(:)*ave_mat(:))
+
+
+ !mask array of mask and region number.
+ DO jj = 1,nreg
+ reg_id_mat(jj) = real(jj-1)
+ mask_id_mat(jj) = real(maskno)
+ END DO
+
+
+ !write text and binary, and note region statistics for current mask for later iom_put
+ IF( lwp ) THEN
+
+ !Write out ascii and binary if requred
+ IF ( ln_diaregmean_bin ) THEN
+ !Writing out regional mean time series to binary files
+ WRITE(numdct_reg_bin) tmp_name,kt,maskno,n_regions_output
+ WRITE(numdct_reg_bin) ave_mat
+ WRITE(numdct_reg_bin) tot_mat
+ WRITE(numdct_reg_bin) var_mat
+ WRITE(numdct_reg_bin) ssq_mat
+ WRITE(numdct_reg_bin) cnt_mat
+ WRITE(numdct_reg_bin) min_mat
+ WRITE(numdct_reg_bin) max_mat
+ ENDIF
+
+ IF ( ln_diaregmean_ascii ) THEN
+ !Writing out regional mean time series to text files
+
+ WRITE(nreg_string, "(I5)") nreg
+ FormatString = "(A30,"//trim(nreg_string)//"F25.3)"
+ WRITE(numdct_reg_txt, FMT="(A30,I6,I6)") tmp_name,kt,maskno
+ WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"ave_mat:", ave_mat
+ WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"tot_mat:", tot_mat
+ WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"var_mat:", var_mat
+ WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"ssq_mat:", ssq_mat
+ WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"cnt_mat:", cnt_mat
+ WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"min_mat:", min_mat
+ WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"max_mat:", max_mat
+ WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"reg_mat:", reg_id_mat
+ WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"msk_mat:", mask_id_mat
+
+ ENDIF
+
+ DO jm = 1,nreg
+ zrmet_ave( reg_ind_cnt) = ave_mat(jm)
+ zrmet_tot( reg_ind_cnt) = tot_mat(jm)
+ zrmet_var( reg_ind_cnt) = var_mat(jm)
+ zrmet_cnt( reg_ind_cnt) = cnt_mat(jm)
+ zrmet_min( reg_ind_cnt) = min_mat(jm)
+ zrmet_max( reg_ind_cnt) = max_mat(jm)
+ zrmet_reg_id( reg_ind_cnt) = reg_id_mat(jm)
+ zrmet_mask_id(reg_ind_cnt) = mask_id_mat(jm)
+
+ reg_ind_cnt = reg_ind_cnt + 1
+ END DO
+
+ ENDIF
+
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean about to deallocated arrays for ',kt,maskno
+ DEALLOCATE(ave_mat,tot_mat,num_mat,var_mat,ssq_mat,cnt_mat,min_mat,max_mat,reg_id_mat,mask_id_mat)
+
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean deallocated arrays for ',kt,maskno
+ IF(lwp)CALL FLUSH(numdct_reg_txt)
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean flushed region mean text for ',kt,maskno
+ END DO
+
+ IF(lwp .AND. verbose) THEN ! Control print
+ WRITE(numout,*) 'dia_regmean ready to start iom_put'
+ CALL FLUSH(numout)
+ ENDIF
+
+ !With current field_def.xml and iodef.xml, these fields must be output, so set to dummy values if not required.
+
+ IF ( ln_diaregmean_nc ) THEN
+
+ zrmet_out(:,:,:) = 0
+ zrmet_val = 0
+ tmp_name_iom = ''
+
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean ready to start iom_put: ',trim(tmp_name)
+
+
+ DO jm = 1,n_regions_output
+ zrmet_val = zrmet_ave(jm)
+! if (zrmet_val .LT. -1e16) zrmet_val = -1e16
+! if (zrmet_val .GT. 1e16) zrmet_val = 1e16
+ if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
+ zrmet_out(:,:,jm) = zrmet_val
+ END DO
+ tmp_name_iom = trim(trim("reg_") // trim(tmp_name) // trim('_ave'))
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)
+ CALL iom_put(trim(tmp_name_iom), zrmet_out )
+ zrmet_out(:,:,:) = 0
+ zrmet_val = 0
+ tmp_name_iom = ''
+
+ DO jm = 1,n_regions_output
+ zrmet_val = zrmet_tot(jm)
+! if (zrmet_val .LT. -1e16) zrmet_val = -1e16
+! if (zrmet_val .GT. 1e16) zrmet_val = 1e16
+ if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
+ zrmet_out(:,:,jm) = zrmet_val
+ END DO
+ tmp_name_iom = trim(trim("reg_") // trim(tmp_name) // trim('_tot'))
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)
+ CALL iom_put( trim(tmp_name_iom), zrmet_out )
+ zrmet_out(:,:,:) = 0
+ zrmet_val = 0
+ tmp_name_iom = ''
+
+ DO jm = 1,n_regions_output
+ zrmet_val = zrmet_var(jm)
+! if (zrmet_val .LT. -1e16) zrmet_val = -1e16
+! if (zrmet_val .GT. 1e16) zrmet_val = 1e16
+ if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
+ zrmet_out(:,:,jm) = zrmet_val
+ END DO
+ tmp_name_iom = trim(trim("reg_") // trim(tmp_name) // trim('_var'))
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)
+ CALL iom_put( trim(tmp_name_iom), zrmet_out )
+ zrmet_out(:,:,:) = 0
+ zrmet_val = 0
+ tmp_name_iom = ''
+
+ DO jm = 1,n_regions_output
+ zrmet_val = zrmet_cnt(jm)
+! if (zrmet_val .LT. -1e16) zrmet_val = -1e16
+! if (zrmet_val .GT. 1e16) zrmet_val = 1e16
+ if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
+ zrmet_out(:,:,jm) = zrmet_val
+ END DO
+ tmp_name_iom = trim(trim("reg_") // trim(tmp_name) // trim('_cnt'))
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)
+ CALL iom_put( trim(tmp_name_iom), zrmet_out )
+ zrmet_out(:,:,:) = 0
+ zrmet_val = 0
+ tmp_name_iom = ''
+
+ DO jm = 1,n_regions_output
+ zrmet_val = zrmet_reg_id(jm)
+! if (zrmet_val .LT. -1e16) zrmet_val = -1e16
+! if (zrmet_val .GT. 1e16) zrmet_val = 1e16
+ if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
+ zrmet_out(:,:,jm) = zrmet_val
+ END DO
+ tmp_name_iom = trim(trim("reg_") // trim(tmp_name) // trim('_reg_id'))
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)
+ CALL iom_put( trim(tmp_name_iom), zrmet_out )
+ zrmet_out(:,:,:) = 0
+ zrmet_val = 0
+ tmp_name_iom = ''
+
+ DO jm = 1,n_regions_output
+ zrmet_val = zrmet_mask_id(jm)
+! if (zrmet_val .LT. -1e16) zrmet_val = -1e16
+! if (zrmet_val .GT. 1e16) zrmet_val = 1e16
+ if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
+ zrmet_out(:,:,jm) = zrmet_val
+ END DO
+ tmp_name_iom = trim(trim("reg_") // trim(tmp_name) // trim('_mask_id'))
+ IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)
+ CALL iom_put( trim(tmp_name_iom), zrmet_out )
+ zrmet_out(:,:,:) = 0
+ zrmet_val = 0
+ tmp_name_iom = ''
+ ELSE
+
+ ALLOCATE( dummy_zrmet(jpi,jpj,n_regions_output), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_wri_region_mean: failed to allocate dummy_zrmet array' )
+
+ DO jm = 1,n_regions_output
+ dummy_zrmet(:,:,jm) = real(jm,wp)
+ END DO
+
+ DO jm = 1,9
+ CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_ave')), dummy_zrmet )
+ CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_tot')), dummy_zrmet )
+ CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_var')), dummy_zrmet )
+ CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_cnt')), dummy_zrmet )
+ CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_reg_id')), dummy_zrmet )
+ CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_mask_id')), dummy_zrmet )
+ END DO
+
+ DEALLOCATE( dummy_zrmet)
+ ENDIF
+
+ DEALLOCATE(zrmet_ave,zrmet_tot,zrmet_var,zrmet_cnt,zrmet_mask_id,zrmet_reg_id,zrmet_min,zrmet_max,zrmet_out)
+
+ IF(lwp .AND. verbose) THEN ! Control print
+ WRITE(numout,*)
+ WRITE(numout,*) 'dia_wri_region_mean finished for ', trim(tmp_name)
+ WRITE(numout,*)
+ CALL FLUSH(numout)
+ ENDIF
+
+ END SUBROUTINE dia_wri_region_mean
+
+
+ !!======================================================================
+END MODULE diaregmean
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90 (revision 10957)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90 (revision 10958)
@@ -11,4 +11,5 @@
USE iom ! I/0 library
USE wrk_nemo ! working arrays
+ USE diaregmean
#if defined key_fabm
USE trc, ONLY: trn
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 (revision 10957)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 (revision 10958)
@@ -47,4 +47,6 @@
USE dia25h ! 25h Mean output
USE diaopfoam ! Diaopfoam output
+ USE diaregmean ! regionalmean
+ USE diapea ! pea
USE iom
USE ioipsl
@@ -405,4 +407,10 @@
IF (ln_diaopfoam) THEN
CALL dia_diaopfoam
+ ENDIF
+ if ( ln_pea ) THEN
+ CALL dia_pea( kt )
+ ENDIF
+ IF (ln_diaregmean) THEN
+ CALL dia_regmean( kt )
ENDIF
!
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90 (revision 10957)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90 (revision 10958)
@@ -106,4 +106,87 @@
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_bnds
!!----------------------------------------------------------------------
+
+
+
+
+
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tmpregion !: temporary region_mask
+ INTEGER, DIMENSION(3) :: zdimsz ! number of elements in each of the 3 dimensions (i.e., lon, lat, no of masks, 297, 375, 4) for an array
+ INTEGER :: zndims ! number of dimensions in an array (i.e. 3, )
+ INTEGER :: inum, nmasks,ierr,maskno,idmaskvar,tmpint
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: tmp_region_mask_real ! tempory region_mask of reals
+
+ LOGICAL :: ln_diaregmean ! region mean calculation
+
+
+ INTEGER :: ios ! Local integer output status for namelist read
+ LOGICAL :: ln_diaregmean_ascii ! region mean calculation ascii output
+ LOGICAL :: ln_diaregmean_bin ! region mean calculation binary output
+ LOGICAL :: ln_diaregmean_nc ! region mean calculation netcdf output
+ LOGICAL :: ln_diaregmean_karamld ! region mean calculation including kara mld terms
+ LOGICAL :: ln_diaregmean_pea ! region mean calculation including pea terms
+ LOGICAL :: ln_diaregmean_diaar5 ! region mean calculation including AR5 SLR terms
+ LOGICAL :: ln_diaregmean_diasbc ! region mean calculation including Surface BC
+
+#if defined key_fabm
+ LOGICAL :: ln_diaregmean_bgc ! region mean calculation including BGC
+#endif
+ ! Read the number region mask to work out how many regions are needed.
+
+#if defined key_fabm
+ NAMELIST/nam_diaregmean/ ln_diaregmean,ln_diaregmean_ascii,ln_diaregmean_bin,ln_diaregmean_nc,&
+ & ln_diaregmean_karamld, ln_diaregmean_pea,ln_diaregmean_diaar5,ln_diaregmean_diasbc,ln_diaregmean_bgc
+#else
+ NAMELIST/nam_diaregmean/ ln_diaregmean,ln_diaregmean_ascii,ln_diaregmean_bin,ln_diaregmean_nc,&
+ & ln_diaregmean_karamld, ln_diaregmean_pea,ln_diaregmean_diaar5,ln_diaregmean_diasbc
+#endif
+
+ ! read in Namelist.
+ !!----------------------------------------------------------------------
+ !
+ REWIND ( numnam_ref ) ! Read Namelist nam_diatmb in referdiatmbence namelist : TMB diagnostics
+ READ ( numnam_ref, nam_diaregmean, IOSTAT=ios, ERR= 901 )
+901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaregmean in reference namelist', lwp )
+
+ REWIND( numnam_cfg ) ! Namelist nam_diatmb in configuration namelist TMB diagnostics
+ READ ( numnam_cfg, nam_diaregmean, IOSTAT = ios, ERR = 902 )
+902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaregmean in configuration namelist', lwp )
+ IF(lwm) WRITE ( numond, nam_diaregmean )
+
+ IF (ln_diaregmean) THEN
+
+ ! Open region mask for region means, and retrieve the size of the mask (number of levels)
+ CALL iom_open ( 'region_mask.nc', inum )
+ idmaskvar = iom_varid( inum, 'mask', kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE.)
+ nmasks = zdimsz(3)
+
+ ! read in the region mask (which contains floating point numbers) into a temporary array of reals.
+ ALLOCATE( tmp_region_mask_real(jpi,jpj,nmasks), STAT= ierr )
+ IF( ierr /= 0 ) CALL ctl_stop( 'dia_regmean_init: failed to allocate tmp_region_mask_real array' )
+
+ ! Use jpdom_unknown to read in a n layer mask.
+ tmp_region_mask_real(:,:,:) = 0
+ CALL iom_get( inum, jpdom_unknown, 'mask', tmp_region_mask_real(1:nlci,1:nlcj,1:nmasks), &
+ & kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nmasks /) )
+
+ CALL iom_close( inum )
+ !Convert the region mask of reals into one of integers.
+
+
+ n_regions_output = 0
+ DO maskno = 1,nmasks
+ tmpint = maxval(int(tmp_region_mask_real(:,:,maskno)))
+ CALL mpp_max( tmpint )
+ n_regions_output = n_regions_output + (tmpint + 1)
+ END DO
+
+
+
+ ELSE
+ n_regions_output = 1
+ ENDIF
+
+
+
#if ! defined key_xios2
ALLOCATE( z_bnds(jpk,2) )
@@ -227,4 +310,11 @@
CALL iom_set_axis_attr( "iax_20C", (/ REAL(20,wp) /) )
CALL iom_set_axis_attr( "iax_28C", (/ REAL(28,wp) /) )
+
+
+
+ CALL iom_set_axis_attr( "region", (/ (REAL(ji,wp), ji=1,n_regions_output) /) )
+
+ CALL iom_set_axis_attr( "noos", (/ (REAL(ji,wp), ji=1,3) /) )
+
! automatic definitions of some of the xml attributs