verticalModes_mod sourceΒΆ

  1module verticalModes_mod
  2  ! MODULE verticalModes_mod (prefix='vms' category='4. Data Object transformations')
  3  !
  4  !:Purpose:  To 1) compute empirical orthogonal functions (EOFs) from either ensemble-derived
  5  !           vertical background-error covariances matrices or a prescribed vertical correlation
  6  !           function (i.e., the so-called vertical modes) and to 2) project back or forth
  7  !           ensemble pertubations onto these modes.
  8  !           Therefore, capablity #2 behaves like a spectral transform but in the vertical
  9  !           dimension.
 10  !
 11  use earthConstants_mod
 12  use utilities_mod
 13  use varNameList_mod
 14  use horizontalCoord_mod
 15  use verticalCoord_mod
 16  use localizationFunction_mod
 17  use calcHeightAndPressure_mod
 18  use gridStatevector_mod
 19  use ensembleStatevector_mod
 20  use timeCoord_mod
 21  use midasMpi_mod
 22  implicit none
 23  save
 24  private
 25  
 26  ! Public derived type
 27  public :: struct_vms
 28  ! Public subroutines
 29  public :: vms_computeModesFromEns, vms_computeModesFromFunction
 30  public :: vms_writeModes, vms_transform
 31
 32  type :: struct_oneVar3d
 33    integer :: nLev
 34    character(len=4) :: varName
 35    real(8), allocatable :: autoCovariance(:,:)
 36    real(8), allocatable :: eigenVectors(:,:)
 37    real(8), allocatable :: eigenValues(:)
 38    real(8), allocatable :: eigenVectorsInv(:,:)
 39  end type struct_oneVar3d
 40  
 41  type :: struct_vms
 42    private
 43    integer :: nVar3d
 44    type(struct_oneVar3d), allocatable :: allVar3d(:)
 45    logical :: initialized = .false.
 46  end type struct_vms
 47
 48contains
 49
 50  !--------------------------------------------------------------------------
 51  ! vms_setup
 52  !--------------------------------------------------------------------------
 53  subroutine vms_setup(varNamesList, vco, vModes)
 54    !
 55    !:Purpose: To setup the vModes structure
 56    !
 57    implicit none
 58
 59    ! Arguments:
 60    character(len=4),          intent(in)    :: varNamesList(:)
 61    type(struct_vco), pointer, intent(in)    :: vco
 62    type(struct_vms),          intent(inout) :: vModes
 63
 64    ! Locals:
 65    integer :: nLev, varNameIndex, var3dIndex
 66    
 67    if (vModes%initialized) return
 68    
 69    !
 70    !- Set parameters 
 71    !
 72    vModes%nVar3d = 0
 73    do varNameIndex = 1, size(varNamesList)
 74      if (vnl_varLevelFromVarname(varNamesList(varNameIndex)) /= 'SF' ) then
 75        vModes%nVar3d = vModes%nVar3d + 1
 76      end if
 77    end do
 78
 79    allocate(vModes%allVar3d(vModes%nVar3d))
 80
 81    var3dIndex = 0
 82    do varNameIndex = 1, size(varNamesList)
 83      if (vnl_varLevelFromVarname(varNamesList(varNameIndex)) /= 'SF' ) then
 84        var3dIndex = var3dIndex + 1
 85        vModes%allVar3d(var3dIndex)%varName = varNamesList(varNameIndex)
 86        if (vnl_varLevelFromVarName(vModes%allVar3d(var3dIndex)%varName) == 'MM') then
 87          nLev = vco_getNumLev(vco,'MM')
 88        else
 89          nLev = vco_getNumLev(vco,'TH')
 90        end if
 91        allocate(vModes%allVar3d(var3dIndex)%autoCovariance(nLev,nLev))
 92        allocate(vModes%allVar3d(var3dIndex)%eigenVectors(nLev,nLev))
 93        allocate(vModes%allVar3d(var3dIndex)%eigenValues(nLev))
 94        allocate(vModes%allVar3d(var3dIndex)%eigenVectorsInv(nLev,nLev))
 95        vModes%allVar3d(var3dIndex)%nLev = nLev
 96      end if
 97    end do
 98
 99    !
100    !- Ending
101    !
102    vModes%initialized = .true.
103
104  end subroutine vms_setup
105  
106  !--------------------------------------------------------------------------
107  ! vms_computeModesFromEns
108  !--------------------------------------------------------------------------
109  subroutine vms_computeModesFromEns(ensPerts,vModes)
110    !
111    !:Purpose: To compute vertical modes from ensemble-derived
112    !           vertical background-error covariances matrices
113    !
114    implicit none
115
116    ! Arguments:
117    type(struct_ens), intent(inout) :: ensPerts
118    type(struct_vms), intent(inout) :: vModes
119
120    ! Locals:
121    type(struct_hco), pointer :: hco_ens
122    type(struct_vco), pointer :: vco_ens
123    real(4),          pointer :: ptr4d_1_r4(:,:,:,:)
124    real(4),          pointer :: ptr4d_2_r4(:,:,:,:)
125    real(8),      allocatable :: latWeight(:) ! Weight given to grid point in the statistic computation
126    real(8),      allocatable :: sumWeight(:,:)
127    integer :: varLevIndex1, varLevIndex2, levIndex1, levIndex2
128    integer :: lonIndex, latIndex, nEns, memberIndex
129    integer :: var3dIndex, myLonBeg, myLonEnd, myLatBeg, myLatEnd
130    character(len=4), pointer :: varNamesList(:)
131    
132    !
133    !- Structure initialization
134    !
135    nullify(varNamesList)
136    call ens_varNamesList(varNamesList, ensPerts)
137
138    vco_ens => ens_getVco(ensPerts)
139    
140    call vms_setup(varNamesList, vco_ens, & ! IN
141                   vModes)                  ! OUT
142
143    !
144    !- Compute the vertical auto-covariance matrices
145    !
146    nEns = ens_getNumMembers(ensPerts)
147    
148    hco_ens => ens_getHco(ensPerts)
149    allocate(latWeight(hco_ens%nj))
150    do latIndex = 1, hco_ens%nj
151      latWeight(latIndex) = cos(hco_ens%lat(latIndex))
152      if (mmpi_myid == 0) then
153        write(*,*) latIndex, hco_ens%lat(latIndex), cos(hco_ens%lat(latIndex))
154      end if
155    end do 
156    
157    call ens_getLatLonBounds(ensPerts, myLonBeg, myLonEnd, myLatBeg, myLatEnd)
158    
159    do var3dIndex = 1, vModes%nVar3d
160
161      allocate(sumWeight(vModes%allVar3d(var3dIndex)%nLev,vModes%allVar3d(var3dIndex)%nLev))
162      sumWeight(:,:)      = 0.d0
163      
164      vModes%allVar3d(var3dIndex)%autoCovariance(:,:) = 0.d0
165
166      do memberIndex = 1, nEns
167
168        !$OMP PARALLEL DO PRIVATE (varLevIndex1,varLevIndex2,levIndex2,levIndex1,ptr4d_1_r4,ptr4d_2_r4,latIndex,lonIndex)
169        do levIndex2 = 1, vModes%allVar3d(var3dIndex)%nLev
170          varLevIndex2 = levIndex2 + ens_getOffsetFromVarName(ensPerts,vModes%allVar3d(var3dIndex)%varName)
171          do levIndex1 = 1, vModes%allVar3d(var3dIndex)%nLev
172            varLevIndex1 = levIndex1 + ens_getOffsetFromVarName(ensPerts,vModes%allVar3d(var3dIndex)%varName)
173
174            ptr4d_1_r4 => ens_getOneLev_r4(ensPerts,varLevIndex1)
175            ptr4d_2_r4 => ens_getOneLev_r4(ensPerts,varLevIndex2)
176          
177            do latIndex = myLatBeg, myLatEnd
178              do lonIndex = myLonBeg, myLonEnd
179
180                vModes%allVar3d(var3dIndex)%autoCovariance(levIndex1,levIndex2) =                     &
181                                     vModes%allVar3d(var3dIndex)%autoCovariance(levIndex1,levIndex2)  &
182                                     + real(ptr4d_1_r4(memberIndex,1,lonIndex,latIndex),8)    &
183                                     * real(ptr4d_2_r4(memberIndex,1,lonIndex,latIndex),8)    &
184                                     * latWeight(latIndex)
185
186                sumWeight(levIndex1,levIndex2) = sumWeight(levIndex1,levIndex2) + latWeight(latIndex)
187                
188              end do
189            end do
190            
191          end do
192        end do
193        !$OMP END PARALLEL DO
194
195      end do ! Loop on Ensemble
196
197      !- Communication
198      call mmpi_allreduce_sumR8_2d(vModes%allVar3d(var3dIndex)%autoCovariance, "GRID")
199      call mmpi_allreduce_sumR8_2d(sumWeight, "GRID")
200
201      !- Apply the appropriate scaling
202      do levIndex2 = 1, vModes%allVar3d(var3dIndex)%nLev
203        do levIndex1 = 1, vModes%allVar3d(var3dIndex)%nLev
204          vModes%allVar3d(var3dIndex)%autoCovariance(levIndex1,levIndex2) =      &
205               vModes%allVar3d(var3dIndex)%autoCovariance(levIndex1,levIndex2)   &
206               / (sumWeight(levIndex1,levIndex2) / dble(nEns) * dble(nEns-1))
207        end do
208      end do
209
210      deallocate(sumWeight)
211
212    end do ! Loop on 3d variables
213
214    !
215    !- Compute the eigenvectors and eigenvalue
216    !
217    call vms_computeModes(vModes)
218
219  end subroutine vms_computeModesFromEns
220
221  !--------------------------------------------------------------------------
222  ! vms_computeModesFromFunction
223  !--------------------------------------------------------------------------
224  subroutine vms_computeModesFromFunction(vco, lengthScaleTop, lengthScaleBot, &
225                                          vModes)
226    !
227    !:Purpose: To compute vertical modes from a prescribed correlation function
228    !
229    implicit none
230
231    ! Arguments:
232    type(struct_vco), pointer, intent(in)    :: vco
233    real(8)         ,          intent(in)    :: lengthScaleTop
234    real(8)         ,          intent(in)    :: lengthScaleBot
235    type(struct_vms),          intent(inout) :: vModes
236
237    ! Locals:
238    real(8), pointer :: vertLocation_MM(:)
239    real(8), pointer :: vertLocation_TH(:)
240    real(8), pointer :: vertLocation(:)
241    real(8) :: pSurfRef, distance, correl
242    real(8) :: lengthScaleGradient
243    real(8) :: lengthScaleLev1, lengthScaleLev2, lengthScaleAvg
244    integer :: levIndex, levIndex1, levIndex2
245    integer :: var3dIndex
246    character(len=4) :: varNamesList(2)
247
248    !
249    !- Structure initialization
250    !
251    varNamesList = (/'P_M','P_T'/)
252    call vms_setup(varNamesList, vco, & ! IN
253                   vModes)              ! OUT
254
255    !
256    !- Compute the vertical auto-covariance matrices
257    !
258    pSurfRef = 101000.D0
259    call czp_fetch1DLevels(vco, pSurfRef,            & ! IN
260                           profM_opt=vertLocation_MM)  ! OUT
261    call czp_fetch1DLevels(vco, pSurfRef,            & ! IN
262                           profT_opt=vertLocation_TH)  ! OUT
263      
264    do levIndex = 1, vco%nLev_M
265      vertLocation_MM(levIndex) = log(vertLocation_MM(levIndex))
266    end do
267    do levIndex = 1, vco%nLev_T
268      vertLocation_TH(levIndex) = log(vertLocation_TH(levIndex))
269    end do
270
271    call lfn_setup('FifthOrder') ! IN
272
273    do var3dIndex = 1, vModes%nVar3d
274
275      if (vnl_varLevelFromVarName(vModes%allVar3d(var3dIndex)%varName) == 'MM') then
276        vertLocation => vertLocation_MM
277      else
278        vertLocation => vertLocation_TH
279      end if
280
281      lengthScaleGradient = (lengthScaleTop-lengthScaleBot) / &
282                            (vertLocation(1)-vertLocation(vModes%allVar3d(var3dIndex)%nLev))
283      
284      do levIndex1 = 1, vModes%allVar3d(var3dIndex)%nLev
285        do levIndex2 = 1, vModes%allVar3d(var3dIndex)%nLev
286          distance = abs(vertLocation(levIndex2) - vertLocation(levIndex1))
287
288          lengthScaleLev1 = lengthScaleGradient * &
289               (vertLocation(levIndex1)-vertLocation(vModes%allVar3d(var3dIndex)%nLev)) + &
290               lengthScaleBot
291          lengthScaleLev2 = lengthScaleGradient * &
292               (vertLocation(levIndex2)-vertLocation(vModes%allVar3d(var3dIndex)%nLev)) + &
293               lengthScaleBot
294          lengthScaleAvg = (lengthScaleLev1+lengthScaleLev2)/2.d0
295
296          if (levIndex1 == 1 .and. mmpi_myid == 0) then
297            write(*,*) 'vms_computeModesFromFunction: function length scale (', &
298                        levIndex2,') = ', lengthScaleAvg
299          end if
300 
301          correl = lfn_response(distance,lengthScaleAvg)
302          vModes%allVar3d(var3dIndex)%autoCovariance(levIndex1,levIndex2) = correl
303        end do
304      end do
305
306    end do
307
308    !
309    !- Compute the eigenvectors and eigenvalue
310    !
311    call vms_computeModes(vModes)
312
313  end subroutine vms_computeModesFromFunction
314
315  !--------------------------------------------------------------------------
316  ! vms_computeModes
317  !--------------------------------------------------------------------------
318  subroutine vms_computeModes(vModes)
319    !
320    !:Purpose: Compute vertical modes from the vertical covariances matrices
321    !           contained in the vModes structure.
322    !
323    implicit none
324
325    ! Arguments:
326    type(struct_vms), intent(inout) :: vModes
327
328    ! Locals:
329    real(8) :: tolerance    
330    integer :: levIndex1, levIndex2, var3dIndex
331    integer :: matrixRank
332    
333    !
334    !- Compute the eigenvectors and eigenvalue
335    !
336    tolerance = 1.0D-50
337    do var3dIndex = 1, vModes%nVar3d
338
339      !- Compute the eigenvectors and eigenvalue
340      call utl_eigenDecomp(vModes%allVar3d(var3dIndex)%autoCovariance, & ! IN
341                           vModes%allVar3d(var3dIndex)%eigenValues,    & ! OUT
342                           vModes%allVar3d(var3dIndex)%eigenVectors,   & ! OUT
343                           tolerance,                                  & ! IN
344                           matrixRank)                                   ! OUT
345      if ( matrixRank < vModes%allVar3d(var3dIndex)%nLev ) then
346        write(*,*) 'varName    =', vModes%allVar3d(var3dIndex)%varName
347        write(*,*) 'matrixRank =', matrixRank
348        write(*,*) 'nLev       =', vModes%allVar3d(var3dIndex)%nLev
349        call utl_abort('vms_computeModes: vertical matrix is rank deficient')
350      end if
351
352      !- Compute the inverse of the eigenVector matrix to go from grid point space
353      !  to vertical mode space
354      !    Note: The inverse of an orthogonal matrix is simply its transpose
355      do levIndex2 = 1, vModes%allVar3d(var3dIndex)%nLev
356        do levIndex1 = 1, vModes%allVar3d(var3dIndex)%nLev
357          vModes%allVar3d(var3dIndex)%eigenVectorsInv(levIndex2,levIndex1) = &
358               vModes%allVar3d(var3dIndex)%eigenVectors(levIndex1,levIndex2)
359        end do
360      end do
361
362    end do
363    
364  end subroutine vms_computeModes
365
366  !--------------------------------------------------------------------------
367  ! vms_transform
368  !--------------------------------------------------------------------------
369  subroutine vms_transform(vModes, VertModesState, GridState, &
370                           TransformDirection, lonBeg, lonEnd, &
371                           latBeg, latEnd, nLev, varName)
372    !
373    !:Purpose: To project back or forth ensemble pertubations onto the vertical
374    !           modes contained in the vModes structure.
375    !
376    implicit none
377
378    ! Arguments:
379    type(struct_vms), intent(in)    :: vModes
380    integer,          intent(in)    :: lonBeg
381    integer,          intent(in)    :: lonEnd
382    integer,          intent(in)    :: latBeg
383    integer,          intent(in)    :: latEnd
384    integer,          intent(in)    :: nLev
385    real(8),          intent(inout) :: VertModesState(lonBeg:lonEnd,latBeg:latEnd,nLev) ! 3D vertical modes coefficients
386    real(8),          intent(inout) :: GridState(lonBeg:lonEnd,latBeg:latEnd,nLev) ! 3D field in grid point space
387    character(len=*), intent(in)    :: TransformDirection ! VertModesToGridPoint or GridPointToVertModes
388    character(len=*), intent(in)    :: varName
389
390    ! Locals:
391    integer :: latIndex, lonIndex, nMode
392    integer :: varIndexAssociated, varIndex
393
394    !
395    !- 1.  Find the appropriate modes for the provided varName
396    !
397    varIndexAssociated = -1
398
399    ! Search for a match in varName
400    do varIndex = 1, vModes%nVar3d
401      if (vModes%allVar3d(varIndex)%varName == varName) then
402        varIndexAssociated = varIndex
403        exit
404      end if
405    end do
406
407    ! Search for a match in level type (MM or TH)
408    if (varIndexAssociated == -1) then
409      if (vnl_varLevelFromVarName(varName) == 'MM') then
410        do varIndex = 1, vModes%nVar3d
411          if (vModes%allVar3d(varIndex)%varName == 'P_M') then
412            varIndexAssociated = varIndex
413            exit
414          end if
415        end do
416      else if (vnl_varLevelFromVarName(varName) == 'TH') then
417        do varIndex = 1, vModes%nVar3d
418          if (vModes%allVar3d(varIndex)%varName == 'P_T') then
419            varIndexAssociated = varIndex
420            exit
421          end if
422        end do
423      else
424        write(*,*) trim(vModes%allVar3d(varIndex)%varName), vnl_varLevelFromVarName(varName)
425        call utl_abort('vms_transform: varName is not on momentum or thermodynamic levels')
426      end if
427    end if
428
429    if (varIndexAssociated == -1) then
430      write(*,*)
431      write(*,*) 'varName in input', trim(varName)
432      write(*,*) 'varNames found in vModes'
433      do varIndex = 1, vModes%nVar3d
434        write(*,*) ' ... ', trim(vModes%allVar3d(varIndex)%varName)
435      end do
436      call utl_abort('vms_transform: could not match varName with any modes')
437    end if
438
439    if (nLev /= vModes%allVar3d(varIndexAssociated)%nLev) then
440      call utl_abort('vms_transform: the number of levels are not consistent')
441    end if
442    nMode = nLev
443    
444    !
445    !- 2.  Do the transform
446    !
447    select case (trim(TransformDirection))
448    case ('GridPointToVertModes')
449
450      !$OMP PARALLEL DO PRIVATE (latIndex,lonIndex)
451      do latIndex = latBeg, latEnd
452        do lonIndex = lonBeg, lonEnd
453          vertModesState(lonIndex,latIndex,1:nMode) = &
454               matmul(vModes%allVar3d(varIndexAssociated)%eigenVectorsInv(1:nMode,1:nLev), &
455                      gridState(lonIndex,latIndex,1:nLev))
456        end do
457      end do
458      !$OMP END PARALLEL DO
459      
460    case ('VertModesToGridPoint')
461
462      !$OMP PARALLEL DO PRIVATE (latIndex,lonIndex)
463      do latIndex = latBeg, latEnd
464        do lonIndex = lonBeg, lonEnd
465          gridState(lonIndex,latIndex,1:nLev) = &
466               matmul(vModes%allVar3d(varIndexAssociated)%eigenVectors(1:nLev,1:nMode), &
467                      vertModesState(lonIndex,latIndex,1:nMode))
468        end do
469      end do
470      !$OMP END PARALLEL DO
471      
472    case default
473      write(*,*)
474      write(*,*) 'Error: TranformDirection Unknown ', trim(TransformDirection)
475      call utl_abort('vms_transform')
476    end select
477    
478  end subroutine vms_transform
479  
480  !--------------------------------------------------------------------------
481  ! vms_writeModes
482  !--------------------------------------------------------------------------
483  subroutine vms_writeModes(vModes)
484    !
485    !:Purpose: To write the content of the provided vModes structure
486    !
487    implicit none
488
489    ! Arguments:
490    type(struct_vms), intent(in) :: vModes
491
492    ! Locals:
493    integer :: var3dIndex, levIndex1, levIndex2
494    integer :: fstouv, fnom, fstfrm, fclos, iunstats, errorID
495    character(len=4), allocatable :: varName3d(:)
496    character(len=128) :: outfilename
497    
498    if (.not. vModes%initialized) then
499      call utl_abort('vms_writeModes: The vModes structure is not initialized')
500    end if
501
502    if (mmpi_myid == 0) then
503      !
504      !- Write in FST file format
505      !
506      allocate(varName3d(vModes%nVar3d))
507      iunstats = 0
508      errorID = fnom(iunstats,'./vertAutoCov.fst','RND',0)
509      errorID = fstouv(iunstats,'RND')      
510      do var3dIndex = 1, vModes%nVar3d
511        call writeMatrix2d_r8(vModes%allVar3d(var3dIndex)%autoCovariance, vModes%allVar3d(var3dIndex)%nlev, &
512                              iunstats, vModes%allVar3d(var3dIndex)%varName,'AUTOCOVAR')
513        call writeMatrix2d_r8(vModes%allVar3d(var3dIndex)%eigenVectors, vModes%allVar3d(var3dIndex)%nlev, &
514                              iunstats, vModes%allVar3d(var3dIndex)%varName,'EIGENVECTORS')
515        call writeArray1d_r8(vModes%allVar3d(var3dIndex)%eigenValues, vModes%allVar3d(var3dIndex)%nlev, &
516                             iunstats, vModes%allVar3d(var3dIndex)%varName,'EIGENVALUES')
517        varName3d(var3dIndex) = vModes%allVar3d(var3dIndex)%varName
518      end do
519      call writeArray1d_c4(varName3d, vModes%nVar3d, iunstats, 'VAR', 'VARNAME')
520      errorID = fstfrm(iunstats)
521      errorID = fclos (iunstats)
522
523      !
524      !- Write in txt format (for plotting purposes)
525      !
526      do var3dIndex = 1, vModes%nVar3d
527        outfilename = "./vModes_autoCovar_"//trim(vModes%allVar3d(var3dIndex)%varName)//".txt"
528        open (unit=99,file=outfilename,action="write",status="new")
529        do levIndex2 = 1, vModes%allVar3d(var3dIndex)%nLev
530          do levIndex1 = 1, vModes%allVar3d(var3dIndex)%nLev
531            if (levIndex1 == vModes%allVar3d(var3dIndex)%nLev) then
532              write(99,'(2X,E10.4)') vModes%allVar3d(var3dIndex)%autoCovariance(levIndex1,levIndex2)
533            else
534              write(99,'(2X,E10.4,$)') vModes%allVar3d(var3dIndex)%autoCovariance(levIndex1,levIndex2)
535            end if
536          end do
537        end do
538        close(unit=99)
539        outfilename = "./vModes_eigenVectors_"//trim(vModes%allVar3d(var3dIndex)%varName)//".txt"
540        open (unit=99,file=outfilename,action="write",status="new")
541        do levIndex2 = 1, vModes%allVar3d(var3dIndex)%nLev
542          do levIndex1 = 1, vModes%allVar3d(var3dIndex)%nLev
543            if (levIndex1 == vModes%allVar3d(var3dIndex)%nLev) then
544              write(99,'(2X,E10.4)') vModes%allVar3d(var3dIndex)%eigenVectors(levIndex1,levIndex2)
545            else
546              write(99,'(2X,E10.4,$)') vModes%allVar3d(var3dIndex)%eigenVectors(levIndex1,levIndex2)
547            end if
548          end do
549        end do
550        close(unit=99)
551        outfilename = "./vModes_eigenValues_"//trim(vModes%allVar3d(var3dIndex)%varName)//".txt"
552        open (unit=99,file=outfilename,action="write",status="new")
553        do levIndex1 = 1, vModes%allVar3d(var3dIndex)%nLev
554          write(99,'(2X,E10.4)') vModes%allVar3d(var3dIndex)%eigenValues(levIndex1)
555        end do
556        close(unit=99)
557      end do
558    end if
559
560  end subroutine vms_writeModes
561
562  !--------------------------------------------------------------------------
563  ! writeMatrix2d_r8
564  !--------------------------------------------------------------------------
565  subroutine writeMatrix2d_r8(matrix2d,rank,iun,nomvar_in,etiket_in)
566    !
567    !:Purpose: To write a 2D matrix rank x rank
568    !
569    implicit none
570
571    ! Arguments:
572    integer,          intent(in) :: rank
573    real(8),          intent(in) :: matrix2d(rank,rank)
574    integer,          intent(in) :: iun
575    character(len=*), intent(in) :: nomvar_in
576    character(len=*), intent(in) :: etiket_in
577
578    ! Locals:
579    real(4), allocatable :: workecr(:,:)
580    real(4) :: work
581    integer :: errorID, fstecr
582    integer :: dateo, npak, ni, nj, nk
583    integer :: ip1, ip2, ip3, deet, npas, datyp
584    integer :: ig1 ,ig2 ,ig3 ,ig4
585    character(len=1)  :: grtyp
586    character(len=4)  :: nomvar
587    character(len=2)  :: typvar
588    character(len=12) :: etiket
589
590    allocate(workecr(rank, rank))
591
592    npak   = -32
593    dateo  = 0
594    deet   = 0
595    npas   = 0
596    ni     = rank
597    nj     = rank
598    nk     = 1
599    ip1    = 0
600    ip2    = 0
601    ip3    = 0
602    typvar = 'XX'
603    nomvar = nomvar_in
604    etiket = etiket_in
605    grtyp  = 'X'
606    ig1    = 0
607    ig2    = 0
608    ig3    = 0
609    ig4    = 0
610    datyp  = 5
611
612    !- Convert to real 4
613    workecr(:,:) = real(matrix2d(:,:),4)
614
615    !- Writing 
616    errorID = fstecr(workecr, work, npak, iun, dateo, deet, npas, ni, nj, &
617                     nk, ip1, ip2, ip3, typvar, nomvar, etiket, grtyp,    &
618                     ig1, ig2, ig3, ig4, datyp, .true.)
619
620    deallocate(workecr)
621
622  end subroutine writeMatrix2d_r8
623
624  !--------------------------------------------------------------------------
625  ! writeArray1d_r8
626  !--------------------------------------------------------------------------
627  subroutine writeArray1d_r8(array,size,iun,nomvar_in,etiket_in)
628    !
629    !:Purpose: To write a 1D array in real 8
630    !
631    implicit none
632
633    ! Arguments:
634    integer,          intent(in) :: size
635    real(8),          intent(in) :: array(size)
636    integer,          intent(in) :: iun
637    character(len=*), intent(in) :: nomvar_in
638    character(len=*), intent(in) :: etiket_in
639
640    ! Locals:
641    real(4), allocatable :: workecr(:)
642    real(4) :: work
643    integer :: errorID, fstecr
644    integer :: dateo, npak, ni, nj, nk
645    integer :: ip1, ip2, ip3, deet, npas, datyp
646    integer :: ig1 ,ig2 ,ig3 ,ig4
647    character(len=1)  :: grtyp
648    character(len=4)  :: nomvar
649    character(len=2)  :: typvar
650    character(len=12) :: etiket
651
652    allocate(workecr(size))
653
654    npak   = -32
655    dateo  = 0
656    deet   = 0
657    npas   = 0
658    ni     = size
659    nj     = 1
660    nk     = 1
661    ip1    = 0
662    ip2    = 0
663    ip3    = 0
664    typvar = 'XX'
665    nomvar = nomvar_in
666    etiket = etiket_in
667    grtyp  = 'X'
668    ig1    = 0
669    ig2    = 0
670    ig3    = 0
671    ig4    = 0
672    datyp  = 5
673
674    !- Convert to real 4
675    workecr(:) = real(array(:),4)
676
677    !- Writing 
678    errorID = fstecr(workecr, work, npak, iun, dateo, deet, npas, ni, nj, &
679                     nk, ip1, ip2, ip3, typvar, nomvar, etiket, grtyp,    &
680                     ig1, ig2, ig3, ig4, datyp, .true.)
681
682    deallocate(workecr)
683    
684  end subroutine writeArray1d_r8
685
686  !--------------------------------------------------------------------------
687  ! writeArray1d_c4
688  !--------------------------------------------------------------------------
689  subroutine writeArray1d_c4(array,size,iun,nomvar_in,etiket_in)
690    !
691    !:Purpose: To write a 1D array in character with len = 4
692    !
693    implicit none
694
695    ! Arguments:
696    integer,          intent(in) :: size
697    character(len=4), intent(in) :: array(size)
698    integer,          intent(in) :: iun
699    character(len=*), intent(in) :: nomvar_in
700    character(len=*), intent(in) :: etiket_in
701
702    ! Locals:
703    real(4) :: work
704    integer :: errorID, fstecr_s
705    integer :: dateo, npak, ni, nj, nk
706    integer :: ip1, ip2, ip3, deet, npas, datyp
707    integer :: ig1 ,ig2 ,ig3 ,ig4
708    character(len=1)  :: grtyp
709    character(len=4)  :: nomvar
710    character(len=2)  :: typvar
711    character(len=12) :: etiket
712
713    npak   = -32
714    dateo  = 0
715    deet   = 0
716    npas   = 0
717    ni     = 4 ! 4 Characters
718    nj     = size
719    nk     = 1
720    ip1    = 0
721    ip2    = 0
722    ip3    = 0
723    typvar = 'XX'
724    nomvar = nomvar_in
725    etiket = etiket_in
726    grtyp  = 'X'
727    ig1    = 0
728    ig2    = 0
729    ig3    = 0
730    ig4    = 0
731    datyp  = 7 ! Character
732
733    !- Writing
734    errorID = fstecr_s(array, work, npak, iun, dateo, deet, npas, ni, nj, &
735                      nk, ip1, ip2, ip3, typvar, nomvar, etiket, grtyp,    &
736                      ig1, ig2, ig3, ig4, datyp, .true.)
737
738  end subroutine writeArray1d_c4
739  
740end module verticalModes_mod