Assumption initial conditions are on z-levels when ln_sco = .true.

[jamesharle]

The initial conditions for runs ln_sco = .true. are assumed to be on z-levels based on gdept_1d. This is not immediately apparent from reading the manual.

A few options would be to either:

• modify code such that depth information is required in the initial condition file. This would allow an initial condition file, comprising a different number of levels and/or gird type, to be read in and interpolated onto gdept_0.
• remove the interpolation loop completely (bearing in mind the on-going simplification process) and impose that the initial conditions read in have to be on gdept_0, requiring a pre-processing step.

[jamesharle]

A fix I used a while back to get me going - could be tidied:

dtatsd.F90

@@ -22 +22 @@
    USE lib_mpp         ! MPP library
    USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
+   USE in_out_manager
+   USE iom
 
    IMPLICIT NONE
    PRIVATE
@@ -33 +35 @@
    LOGICAL , PUBLIC ::   ln_tsd_tradmp    !: internal damping toward input data flag
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read)
+#if defined key_gen_IC
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,: ) ::   gdept_init 
+#endif
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -147 +152 @@
       !
       INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies
       INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1   ! local integers
+#if defined key_gen_IC
+      INTEGER ::   id, jpk_init, jpk_initm1   ! local integers
+      INTEGER , DIMENSION(4)            ::   ddims
+#endif
       REAL(wp)::   zl, zi
       REAL(wp), POINTER, DIMENSION(:) ::  ztp, zsp   ! 1D workspace
       !!----------------------------------------------------------------------
@@ -154 +163 @@
       IF( nn_timing == 1 )  CALL timing_start('dta_tsd')
       !
       CALL fld_read( kt, 1, sf_tsd )      !==   read T & S data at kt time step   ==!
+
+#if defined key_gen_IC
+      CALL iom_open ( sf_tsd(jp_tem)%clname, sf_tsd(jp_tem)%num ) 
+      !! get dimensions
+      id = iom_varid( sf_tsd(jp_tem)%num, 'gdept', ddims )
+      jpk_init = ddims(3)
+      IF(lwp) WRITE(numout,*) 'Dimensions of ICs: ', ddims, jpk_init
+      ALLOCATE( gdept_init(jpi,jpj,jpk_init) )
+      CALL iom_get ( sf_tsd(jp_tem)%num, jpdom_data, 'gdept', gdept_init,1)
+      CALL iom_close( sf_tsd(jp_tem)%num )   ! Close the input file
+      jpk_initm1=jpk_init-1
+#endif
       !
       !
       !                                   !==   ORCA_R2 configuration and T & S damping   ==! 
@@ -221 +242 @@
             DO ji = 1, jpi
                DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
                   zl = gdept_0(ji,jj,jk)
+#if defined key_gen_IC
+                  IF(     zl < gdept_init(ji,jj,1  ) ) THEN          ! above the first level of data
+#else
                   IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data
+#endif
                      ztp(jk) =  ptsd(ji,jj,1    ,jp_tem)
                      zsp(jk) =  ptsd(ji,jj,1    ,jp_sal)
+#if defined key_gen_IC
+                  ELSEIF( zl > gdept_init(ji,jj,jpk_init) ) THEN          ! below the last level of data
+#else
                   ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data
+#endif
                      ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem)
                      zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal)
                   ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
+#if defined key_gen_IC
+                     DO jkk = 1, jpk_initm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
+                        IF( (zl-gdept_init(ji,jj,jkk)) * (zl-gdept_init(ji,jj,jkk+1)) <= 0._wp ) THEN
+                           zi = ( zl - gdept_init(ji,jj,jkk) ) / (gdept_init(ji,jj,jkk+1)-gdept_init(ji,jj,jkk))
+#else
                      DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
                         IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
                            zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
+#endif
                            ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi 
                            zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi
                         ENDIF
@@ -303 +338 @@
                                         DEALLOCATE( sf_tsd(jp_sal)%fnow )     ! S arrays in the structure
          IF( sf_tsd(jp_sal)%ln_tint )   DEALLOCATE( sf_tsd(jp_sal)%fdta )
                                         DEALLOCATE( sf_tsd              )     ! the structure itself
+#if defined key_gen_IC
+         DEALLOCATE( gdept_init )
+#endif
       ENDIF
       !
       IF( nn_timing == 1 )  CALL timing_stop('dta_tsd')

[jamesharle]

I'm currently writing this code into my copy of the trunk. I was also thinking of removing this loop and putting it into domain_cfg TOOLS so the initial conditions are mapped onto the native vertical grid and written out when domain_cfg is used. That way they can be read in a runtime without the need for on-the-fly interpolation.

[mathiot]

James, one other solution could be to clarify the documentation and if you really want to to implement the capability and an action in the proposal for 2020 WP.

[jamesharle]

OK - I have an updated changeset for v4.0 that I've been using for ORCHESTRA, so won't be too much effort to add. I'll look over the documentation and discuss with Andrew.