[2136] | 1 | PROGRAM create_coordinates |
---|
| 2 | !!----------------------------------------------------------- |
---|
| 3 | !! |
---|
| 4 | !! to create regional coordinates file (e.g. 1_coordinates.nc) |
---|
| 5 | !! for fine grid from coarse grid (e.g. coordinates.nc) |
---|
| 6 | !! |
---|
| 7 | !! to make it we use a 4th order polynomial interpolation |
---|
| 8 | !! |
---|
| 9 | !! Created by Brice Lemaire on 12/2009. |
---|
| 10 | !! |
---|
| 11 | !!----------------------------------------------------------- |
---|
| 12 | USE netcdf |
---|
| 13 | USE cfg_tools |
---|
| 14 | USE mixed_grid |
---|
| 15 | USE domain |
---|
| 16 | !------------ |
---|
| 17 | IMPLICIT NONE |
---|
| 18 | ! |
---|
| 19 | INTEGER :: narg,iargc |
---|
| 20 | INTEGER :: nstat, imodulo |
---|
| 21 | INTEGER :: nik, njk, npk |
---|
| 22 | INTEGER :: ival |
---|
| 23 | CHARACTER(len=80) :: cnmlname, cfingrdname |
---|
| 24 | ! |
---|
| 25 | ! |
---|
| 26 | ! |
---|
| 27 | !************************************* |
---|
| 28 | ! Read input file (namelist file) |
---|
| 29 | ! using types.f90 |
---|
| 30 | !************************************* |
---|
| 31 | narg = iargc() |
---|
| 32 | ! |
---|
| 33 | IF (narg == 0) THEN |
---|
| 34 | cnmlname = 'namelist.input' |
---|
| 35 | ELSE |
---|
| 36 | CALL getarg(1,cnmlname) |
---|
| 37 | ENDIF |
---|
| 38 | ! |
---|
| 39 | CALL read_namelist(cnmlname) |
---|
| 40 | ! |
---|
| 41 | nstat = read_coordinates(TRIM(cn_parent_coordinate_file),scoagrd) |
---|
| 42 | IF (nstat/=1) THEN |
---|
| 43 | WRITE(*,*) 'unable to open netcdf file : ',cn_parent_coordinate_file |
---|
| 44 | STOP |
---|
| 45 | ENDIF |
---|
| 46 | ! |
---|
| 47 | nsizex = SIZE(scoagrd%glamt,1) |
---|
| 48 | nsizey = SIZE(scoagrd%glamt,2) |
---|
| 49 | ! |
---|
| 50 | WRITE(*,*) '' |
---|
| 51 | WRITE(*,*) 'Size of input matrix: ' |
---|
| 52 | WRITE(*,*) '(',nsizex,';', nsizey,')' |
---|
| 53 | WRITE(*,*) '' |
---|
| 54 | ! |
---|
| 55 | WRITE(*,*) 'Domain: ' |
---|
| 56 | WRITE(*,*) '' |
---|
| 57 | WRITE(*,*) ' ', ' min(1,1:nsizey) ', ' max(nsizex,1:nsizey) ' |
---|
| 58 | WRITE(*,*) 'longitude: ', MINVAL(scoagrd%glamt(1,1:nsizey)), ' --> ', MAXVAL(scoagrd%glamf(nsizex,1:nsizey)) |
---|
| 59 | WRITE(*,*) 'latitude: ', MINVAL(scoagrd%gphit(1:nsizex,1)), ' --> ', MAXVAL(scoagrd%gphif(1:nsizex,nsizey)) |
---|
| 60 | WRITE(*,*) '' |
---|
| 61 | ! |
---|
| 62 | !************************************* |
---|
| 63 | ! Check the values read in the namelist |
---|
| 64 | !************************************* |
---|
| 65 | IF(nn_imin.LE.0.OR.nn_imin.GT.nsizex)THEN |
---|
| 66 | WRITE(*,*) 'Wrong values in namelist:' |
---|
| 67 | WRITE(*,*) 'nn_imin must be greater than 0 and less than', nsizex |
---|
| 68 | STOP |
---|
| 69 | ENDIF |
---|
| 70 | ! |
---|
| 71 | IF(nn_imax.LE.0.OR.nn_imax.GT.nsizex)THEN |
---|
| 72 | WRITE(*,*) 'Wrong values in namelist:' |
---|
| 73 | WRITE(*,*) 'nn_imax must be greater than 0 and less than', nsizex |
---|
| 74 | STOP |
---|
| 75 | ENDIF |
---|
| 76 | ! |
---|
| 77 | IF(nn_jmin.LE.0.OR.nn_jmin.GT.nsizey)THEN |
---|
| 78 | WRITE(*,*) 'Wrong values in namelist:' |
---|
| 79 | WRITE(*,*) 'nn_jmin must be greater than 0 and less than', nsizey |
---|
| 80 | STOP |
---|
| 81 | ENDIF |
---|
| 82 | ! |
---|
| 83 | IF(nn_jmax.LE.0.OR.nn_jmax.GT.nsizey)THEN |
---|
| 84 | WRITE(*,*) 'Wrong values in namelist:' |
---|
| 85 | WRITE(*,*) 'nn_jmax must be greater than 0 and less than', nsizey |
---|
| 86 | STOP |
---|
| 87 | ENDIF |
---|
| 88 | ! |
---|
| 89 | !************************************* |
---|
| 90 | ! Define the sub-domain chosen by user |
---|
| 91 | !************************************* |
---|
| 92 | WRITE(*,*) 'Domain defined by user: ' |
---|
| 93 | WRITE(*,*) '' |
---|
| 94 | WRITE(*,*) ' ', ' min ', ' max ' |
---|
| 95 | WRITE(*,*) 'longitude: ', MINVAL(scoagrd%glamt(nn_imin,nn_jmin:nn_jmax)), ' --> ', MAXVAL(scoagrd%glamf(nn_imax,nn_jmin:nn_jmax)) |
---|
| 96 | WRITE(*,*) 'latitude: ', MINVAL(scoagrd%gphit(nn_imin:nn_imax,nn_jmin)), ' --> ', MAXVAL(scoagrd%gphif(nn_imin:nn_imax,nn_jmax)) |
---|
| 97 | WRITE(*,*) '' |
---|
| 98 | ! |
---|
| 99 | IF(nn_imin.EQ.nn_imax.AND.nn_jmin.EQ.nn_jmax) THEN |
---|
| 100 | nglobal = .TRUE. |
---|
| 101 | WRITE(*,*) 'Size of domain: GLOBAL' |
---|
| 102 | ELSE |
---|
| 103 | nglobal = .FALSE. |
---|
| 104 | WRITE(*,*) 'Size of domain: REGIONAL' |
---|
| 105 | ENDIF |
---|
| 106 | ! |
---|
| 107 | IF(cn_position_pivot.EQ.'F-grid') THEN !case of ORCA05 & ORCA025 grids |
---|
| 108 | npivot = 0 |
---|
| 109 | ELSEIF(cn_position_pivot.EQ.'T-grid') THEN !case of ORCA2 grid |
---|
| 110 | npivot = 1 |
---|
| 111 | ENDIF |
---|
| 112 | ! |
---|
| 113 | nmid = nsizex/2 + npivot |
---|
| 114 | ! |
---|
| 115 | ! |
---|
| 116 | ! |
---|
| 117 | !************************************* |
---|
| 118 | ! Redefine for particular cases |
---|
| 119 | !************************************* |
---|
| 120 | IF(((nn_imax - nn_imin).EQ.-1).OR.((nn_imax - nn_imin).EQ.-2)) THEN |
---|
| 121 | nn_imax = nn_imin |
---|
| 122 | ELSEIF(ABS(nn_imin - nn_imax).GE.nsizex-1) THEN |
---|
| 123 | nn_imax = nn_imin |
---|
| 124 | ELSEIF((nn_imin.EQ.1).AND.(nn_imax.EQ.nsizex)) THEN |
---|
| 125 | nn_imax = nn_imin |
---|
| 126 | ELSEIF(nn_imin.EQ.1.AND.nn_imax.GT.nn_imin) THEN |
---|
| 127 | nn_imin = nsizey-2 |
---|
| 128 | ELSEIF(nn_imax.EQ.nsizex) THEN |
---|
| 129 | nn_imax = 3 |
---|
| 130 | ELSEIF(nn_imin.EQ.nsizex) THEN |
---|
| 131 | nn_imin = nsizex-1 |
---|
| 132 | ELSEIF(nn_imax.EQ.1.AND.nn_imax.NE.nn_imin) THEN |
---|
| 133 | nn_imax = 2 |
---|
| 134 | ENDIF |
---|
| 135 | ! |
---|
| 136 | ! |
---|
| 137 | ! |
---|
| 138 | !************************************* |
---|
| 139 | ! We want to fix equator along T and U-points |
---|
| 140 | !************************************* |
---|
| 141 | imodulo = MOD(nn_rhoy,2) |
---|
| 142 | ! |
---|
| 143 | IF(.NOT.nglobal.AND.(scoagrd%gphit(nn_imin,nn_jmin)*scoagrd%gphit(nn_imin,nn_jmax)).LT.0.AND.imodulo.EQ.0) THEN |
---|
| 144 | nequator = 1 |
---|
| 145 | ELSEIF(nglobal.AND.imodulo.EQ.0) THEN |
---|
| 146 | nequator = 1 |
---|
| 147 | ELSE |
---|
| 148 | nequator = 0 |
---|
| 149 | ENDIF |
---|
| 150 | ! |
---|
| 151 | ! |
---|
| 152 | ! |
---|
| 153 | !************************************* |
---|
| 154 | ! CREATE MIXED GRID |
---|
| 155 | !************************************* |
---|
| 156 | CALL define_domain |
---|
| 157 | ! |
---|
| 158 | CALL define_mixed_grid |
---|
| 159 | ! |
---|
| 160 | ! |
---|
| 161 | ! |
---|
| 162 | !************************************* |
---|
| 163 | !!! CALCULATE FINE GRID DIMENSIONS |
---|
| 164 | !************************************* |
---|
| 165 | IF(.NOT.nglobal) THEN |
---|
| 166 | IF(nn_rhox.EQ.1) THEN |
---|
| 167 | nxfine = nxcoag |
---|
| 168 | ELSE |
---|
| 169 | nxfine = (nxcoag-2)*nn_rhox + 1 |
---|
| 170 | ENDIF |
---|
| 171 | ! |
---|
| 172 | IF(nn_rhoy.EQ.1) THEN |
---|
| 173 | nyfine = nycoag |
---|
| 174 | ELSE |
---|
| 175 | nyfine = (nycoag-2)*nn_rhoy + 1 |
---|
| 176 | ENDIF |
---|
| 177 | ! |
---|
| 178 | ELSEIF(nglobal) THEN |
---|
| 179 | IF(nn_rhox.GT.1) THEN |
---|
| 180 | nxfine = (nxgmix - 4*(nn_rhox-1))/2 |
---|
| 181 | ELSEIF(nn_rhox.EQ.1) THEN |
---|
| 182 | nxfine = nsizex |
---|
| 183 | ENDIF |
---|
| 184 | ! |
---|
| 185 | IF(nn_rhoy.GT.1) THEN |
---|
| 186 | nyfine = (nygmix - (2*nn_rhoy))/2 - nn_rhoy + nequator |
---|
| 187 | ELSEIF(nn_rhoy.EQ.1) THEN |
---|
| 188 | nyfine = nsizey |
---|
| 189 | ENDIF |
---|
| 190 | ! |
---|
| 191 | ENDIF |
---|
| 192 | ! |
---|
| 193 | WRITE(*,*) '' |
---|
| 194 | WRITE(*,*) '*** SIZE OF FINE GRID ***' |
---|
| 195 | WRITE(*,*) nxfine, ' x ', nyfine |
---|
| 196 | WRITE(*,*) '' |
---|
| 197 | ! |
---|
| 198 | ! |
---|
| 199 | ! |
---|
| 200 | !************************************* |
---|
| 201 | ! Interpolation inside the mixed grid |
---|
| 202 | ! from cfg_tools.f90 |
---|
| 203 | !************************************* |
---|
| 204 | IF(nn_rhox.GT.1.OR.nn_rhoy.GT.1) THEN |
---|
| 205 | CALL interp_grid |
---|
| 206 | ENDIF |
---|
| 207 | ! |
---|
| 208 | ! |
---|
| 209 | ! |
---|
| 210 | !************************************* |
---|
| 211 | ! Define name of child coordinate file |
---|
| 212 | ! coordinates.nc -> 1_coordinates.nc |
---|
| 213 | !************************************* |
---|
| 214 | CALL set_child_name(cn_parent_coordinate_file,cfingrdname) |
---|
| 215 | ! |
---|
| 216 | ! |
---|
| 217 | ! |
---|
| 218 | !************************************* |
---|
| 219 | ! Allocation of child grid elements |
---|
| 220 | ! from types.f90 |
---|
| 221 | !************************************* |
---|
| 222 | CALL grid_allocate(sfingrd,nxfine,nyfine) |
---|
| 223 | ! |
---|
| 224 | ! |
---|
| 225 | ! |
---|
| 226 | !************************************* |
---|
| 227 | ! Break the mixed grid smixgrd into 4 grids |
---|
| 228 | ! from cfg_tools.f90 |
---|
| 229 | !************************************* |
---|
| 230 | CALL child_grid |
---|
| 231 | ! |
---|
| 232 | ! |
---|
| 233 | ! |
---|
| 234 | !************************************* |
---|
| 235 | ! Read parent coordinate file |
---|
| 236 | ! from readwrite.f90 |
---|
| 237 | !************************************* |
---|
| 238 | nstat = write_coordinates(cfingrdname,sfingrd,nxfine,nyfine) |
---|
| 239 | IF (nstat/=1) THEN |
---|
| 240 | WRITE(*,*)"unable to write netcdf file : ",cfingrdname |
---|
| 241 | STOP |
---|
| 242 | ENDIF |
---|
| 243 | ! |
---|
| 244 | CALL grid_deallocate(scoagrd) |
---|
| 245 | CALL grid_deallocate(sfingrd) |
---|
| 246 | CALL mixed_grid_deallocate(smixgrd) |
---|
| 247 | ! |
---|
| 248 | ! |
---|
| 249 | ! |
---|
| 250 | END PROGRAM create_coordinates |
---|