source: trunk/LMDZ.MARS/libf/phymars/surfini.F @ 1541

Last change on this file since 1541 was 1541, checked in by emillour, 9 years ago

Mars GCM:

  • fix for 1D in writediagfi to enable writing at "ecritphy" rate.
  • removed iniprint.h from phymars/dyn1d since it is in "misc"
  • Some code cleanup in anticipation of future updates:
    • changed variable names in comgeomphy.F90: give them more explicit names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
    • removed long(), lati() and area() from comgeomfi_h.F90, use longitude(), latitude() and cell_are() from comgeomphy.F90 instead

EM
File size: 20.2 KB
Line 
1      SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb)
2   ! to use  'getin'
3      USE ioipsl_getincom
4      use netcdf
5      use tracer_mod, only: nqmx, noms
6      use comgeomphy, only: longitude, latitude ! in radians
7      use surfdat_h, only: watercaptag, frost_albedo_threshold,
8     &                     albedo_h2o_ice, inert_h2o_ice, albedodat,
9     &                     albedice, dryness
10#ifndef MESOSCALE
11      use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
12      use mod_phys_lmdz_para, only : is_master, gather, scatter
13#endif
14      USE comcstfi_h
15      use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
16      IMPLICIT NONE
17c=======================================================================
18c
19c   creation des calottes pour l'etat initial
20c
21c=======================================================================
22c-----------------------------------------------------------------------
23c   Declarations:
24c   -------------
25#include "callkeys.h"
26#include "datafile.h"
27
28      integer,intent(in) :: ngrid ! number of atmospheric columns
29      real,intent(in) :: piceco2(ngrid) ! CO2 ice thickness
30      real,intent(inout) :: qsurf(ngrid,nqmx) ! tracer on surface (kg/m2)
31      real,intent(out) :: psolaralb(ngrid,2) ! albedo
32
33      INTEGER ig,icap,iq,alternate
34      REAL icedryness ! ice dryness
35     
36      ! longwatercaptag is watercaptag. Trick for some compilers
37      LOGICAL, DIMENSION(100000) :: longwatercaptag
38     
39! There are 3 different modes for ice distribution:
40! icelocationmode = 1 ---> based on data from surface.nc
41! icelocationmode = 2 ---> directly predefined for GCM resolutions 32x24 or 64x48
42! icelocationmode = 3 ---> based on logical relations for latitude and longitude
43! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
44      INTEGER,SAVE :: icelocationmode = 2
45       
46       
47      !in case icelocationmode == 1
48      INTEGER i,j
49      INTEGER     imd,jmd
50      PARAMETER   (imd=360,jmd=180)
51      REAL        zdata(imd,jmd)
52      REAL        zelat,zelon
53
54#ifndef MESOSCALE
55      INTEGER nb_ice(klon_glo,2)   ! number of counts | detected ice for GCM grid
56#endif
57      INTEGER latice(nbp_lat-1,2),lonice (nbp_lon,2) ! number of counts | detected ice along lat & lon axis
58
59      REAL step,count,ratiolat
60
61      INTEGER   ierr,nid,nvarid
62     
63      REAL,SAVE :: min_icevalue = 500.
64      character(len=50) :: string = 'thermal'
65     
66      character (len=100) :: zedatafile
67
68#ifdef MESOSCALE
69
70      do ig=1,ngrid
71
72         !write(*,*) "all qsurf to zero. dirty."
73         do iq=1,nqmx
74         qsurf(ig,iq)=0.  !! on jette les inputs GCM
75                          !! on regle juste watercaptag
76                          !! il faudrait garder les inputs GCM
77                          !! si elles sont consequentes
78         enddo
79         if ( ( latitude(ig)*180./pi .gt. 70. ) .and.
80     .        ( albedodat(ig) .ge. 0.26   ) )  then
81                 write(*,*)"outlier ",ig
82                 watercaptag(ig)  = .true.
83                 dryness(ig)      = 1.
84                 albedodat(ig)    = albedo_h2o_ice  !! pour output
85         else
86                 watercaptag(ig)  = .false.
87                 dryness(ig)      = 1.
88         endif
89
90      enddo
91#endif
92! problem with nested precompiling flags
93
94#ifndef MESOSCALE
95      ! to handle parallel cases
96#if CPP_PARA
97      logical watercaptag_glo(klon_glo)
98      real dryness_glo(klon_glo)
99      real lati_glo(klon_glo)
100      real long_glo(klon_glo)
101#else
102      logical watercaptag_glo(ngrid)
103      real dryness_glo(ngrid)
104      real lati_glo(ngrid)
105      real long_glo(ngrid)
106#endif
107#endif
108
109#ifndef MESOSCALE
110
111c
112c=======================================================================
113! Initialize watercaptag (default is false)
114      watercaptag_glo(:)=.false.
115
116c     water ice outliers
117c     ------------------------------------------
118
119      IF ((water) .and. (caps)) THEN
120     
121c Perennial H20 north cap defined by watercaptag=true (allows surface to be
122c hollowed by sublimation in vdifc).
123
124c We might not want albedodat to be modified because it is used to write
125c restart files. Instead, albedo is directly modified when needed (i.e.
126c if we have watercaptag and no co2 ice), below and in albedocaps.F90
127
128c       "Dryness coefficient" controlling the evaporation and
129c        sublimation from the ground water ice (close to 1)
130c        HERE, the goal is to correct for the fact
131c        that the simulated permanent water ice polar caps
132c        is larger than the actual cap and the atmospheric
133c        opacity not always realistic.
134
135         alternate = 0
136         
137         if (ngrid .ne. 1) then
138           watercaptag(:) = .false.
139           longwatercaptag(:) = .false.
140         endif
141         
142         write(*,*) "surfini: Ice dryness ?"
143         icedryness=1. ! default value
144         call getin("icedryness",icedryness)
145         write(*,*) "surfini: icedryness = ",icedryness
146         dryness (:) = icedryness
147         
148      ! To be able to run in parallel, we work on the full grid
149      ! and dispatch results afterwards
150
151      ! start by geting latitudes and logitudes on full grid
152      ! (in serial mode, this is just a copy)
153      call gather(latitude,lati_glo)
154      call gather(longitude,long_glo)
155
156      if (is_master) then
157
158        IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
159     
160         print*, 'ngrid = 1, do no put ice caps in surfini.F'
161
162        ELSE IF (icelocationmode .eq. 1) THEN
163     
164         print*,'Surfini: ice caps defined from surface.nc'
165           
166! This method detects ice as gridded value above min_icevalue in the field "string" from surface.nc
167! Typically, it is for thermal inertia above 500 tiu.
168! Two conditions are verified:
169! 1. GCM ice caps are defined such as area is conserved for a given latitude
170! (the approximation is that all points within the GCM latitude resolution have the same area).
171! 2. caps are placed to fill the GCM points with the most detected ice first.
172     
173
174           
175         zedatafile = trim(datafile)
176 
177       
178         ierr=nf90_open(trim(zedatafile)//'/surface.nc',
179     &   NF90_NOWRITE,nid)
180     
181         IF (ierr.NE.nf90_noerr) THEN
182       write(*,*)'Error : cannot open file surface.nc '
183       write(*,*)'(in phymars/surfini.F)'
184       write(*,*)'It should be in :',trim(zedatafile),'/'
185       write(*,*)'1) You can set this path in the callphys.def file:'
186       write(*,*)'   datadir=/path/to/the/datafiles'
187       write(*,*)'2) If necessary, surface.nc (and other datafiles)'
188       write(*,*)'   can be obtained online on:'
189       write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
190       CALL ABORT
191         ENDIF
192     
193     
194         ierr=nf90_inq_varid(nid, string, nvarid)
195         if (ierr.ne.nf90_noerr) then
196          write(*,*) 'surfini error, cannot find ',trim(string)
197          write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
198          write(*,*)trim(nf90_strerror(ierr))
199          stop
200         endif
201
202         ierr=nf90_get_var(nid, nvarid, zdata)
203
204         if (ierr.ne.nf90_noerr) then
205          write(*,*) 'surfini: error failed loading ',trim(string)
206          write(*,*)trim(nf90_strerror(ierr))
207          stop
208         endif
209 
210                     
211         ierr=nf90_close(nid)
212 
213
214         nb_ice(:,1) = 1 ! default: there is no ice
215         latice(:,1) = 1
216         lonice(:,1) = 1
217         nb_ice(:,2) = 0
218         latice(:,2) = 0
219         lonice(:,2) = 0
220         !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
221
222         ! loop over the GCM grid - except for poles (ig=1 and ngrid)
223         do ig=2,klon_glo-1
224     
225          ! loop over the surface file grid     
226          do i=1,imd
227           do j=1,jmd
228             zelon = i - 180.
229             zelat = 90. - j
230            if ((abs(lati_glo(ig)*180./pi-zelat).le.
231     &           90./real(nbp_lat-1)) .and.
232     &          (abs(long_glo(ig)*180./pi-zelon).le.
233     &           180./real(nbp_lon))) then
234              ! count all points in that GCM grid point
235              nb_ice(ig,1) = nb_ice(ig,1) + 1
236              if (zdata(i,j) > min_icevalue)
237                 ! count all detected points in that GCM grid point
238     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
239             endif
240           enddo
241          enddo 
242
243        ! projection of nb_ice on GCM lat and lon axes
244          latice(1+(ig-2)/nbp_lon,:) =
245     &     latice(1+(ig-2)/nbp_lon,:) + nb_ice(ig,:)
246          lonice(1+mod(ig-2,nbp_lon),:) =
247     &     lonice(1+mod(ig-2,nbp_lon),:) + nb_ice(ig,:) ! lonice is USELESS ...
248
249         enddo ! of do ig=2,klon_glo-1
250     
251
252     
253         ! special case for poles
254         nb_ice(1,2)   = 1  ! ice prescribed on north pole
255         latice(1,:)   = nb_ice(1,:)
256         lonice(1,:)   = nb_ice(1,:)
257         latice(nbp_lat-1,:) = nb_ice(ngrid,:)
258         lonice(nbp_lon,:) = nb_ice(ngrid,:)
259     
260     
261!      print*, 'latice TOT', latice(:,1)
262!      print*, 'latice FOUND', latice(:,2)
263!      print*, 'lonice TOT', lonice(:,1)
264!      print*, 'lonice FOUND', lonice(:,2)
265     
266!      print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
267!      print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
268     
269!      print*,''
270!      print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
271!      print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
272     
273   
274         ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
275         do i=1,(nbp_lat-1)/2
276          step  = 1. ! threshold to add ice cap
277          count = 0. ! number of ice GCM caps at this latitude
278          ! ratiolat is the ratio of area covered by ice within this GCM latitude range
279          ratiolat  = real(latice(i,2))/real(latice(i,1))
280          !print*,'i',i,(i-1)*iim+2,i*iim+1
281     
282          ! put ice caps while there is not enough ice,
283          ! as long as the threshold is above 20%
284          do while ((count.le.ratiolat*nbp_lon).and.(step.ge.0.2))
285           count = 0.
286           ! loop over GCM longitudes
287           do j=1,nbp_lon
288            ! if the detected ice ratio in the GCM grid point
289            ! is more than 'step', then add ice
290            if (real(nb_ice((i-1)*nbp_lon+1+j,2))
291     &        / real(nb_ice((i-1)*nbp_lon+1+j,1)) .ge. step) then
292                  watercaptag_glo((i-1)*nbp_lon+1+j) = .true.
293                  count = count + 1
294            endif
295           enddo ! of do j=1,nbp_lon
296           !print*, 'step',step,count,ratiolat*nbp_lon
297           step = step - 0.01
298          enddo ! of do while
299          !print*, 'step',step,count,ratiolat*nbp_lon
300
301         enddo ! of do i=1,jjm/2
302           
303
304        ELSE IF (icelocationmode .eq. 2) THEN
305     
306         print*,'Surfini: predefined ice caps'
307     
308         if ((nbp_lon.eq.32).and.((nbp_lat-1).eq.24)) then ! 32x24
309           
310          print*,'water ice caps distribution for 32x24 resolution'
311          longwatercaptag(1:9)    = .true. ! central cap - core
312          longwatercaptag(26:33)  = .true. ! central cap
313          longwatercaptag(1:33)  = .true. ! central cap
314          longwatercaptag(56)  = .true. ! central cap
315          longwatercaptag(58)  = .true. ! central cap
316          longwatercaptag(60)  = .true. ! central cap
317          longwatercaptag(62)  = .true. ! central cap
318          longwatercaptag(64)  = .true. ! central cap
319!---------------------   OUTLIERS  ----------------------------
320
321         else if ((nbp_lon.eq.64).and.((nbp_lat-1).eq.48)) then ! 64x48
322
323          print*,'water ice caps distribution for 64x48 resolution'
324          longwatercaptag(1:65)   = .true. ! central cap - core
325          longwatercaptag(75:85)  = .true. ! central cap
326          longwatercaptag(93:114) = .true. ! central cap
327!---------------------   OUTLIERS  ----------------------------
328          if (.true.) then
329          longwatercaptag(136)    = .true. ! outlier, lat = 78.75
330          longwatercaptag(138)    = .true. ! outlier, lat = 78.75
331          longwatercaptag(140)    = .true. ! outlier, lat = 78.75
332          longwatercaptag(142)    = .true. ! outlier, lat = 78.75
333          longwatercaptag(161)    = .true. ! outlier, lat = 78.75
334          longwatercaptag(163)    = .true. ! outlier, lat = 78.75
335          longwatercaptag(165)    = .true. ! outlier, lat = 78.75
336          longwatercaptag(183)    = .true. ! outlier, lat = 78.75
337          longwatercaptag(185)    = .true. ! outlier, lat = 78.75
338          longwatercaptag(187)    = .true. ! outlier, lat = 78.75
339          longwatercaptag(189)    = .true. ! outlier, lat = 78.75
340          longwatercaptag(191)    = .true. ! outlier, lat = 78.75
341          longwatercaptag(193)    = .true. ! outlier, lat = 78.75
342          longwatercaptag(194)    = .true. ! outlier, lat = 75
343          longwatercaptag(203)    = .true. ! outlier, lat = 75
344          longwatercaptag(207)    = .true. ! outlier, lat = 75
345          longwatercaptag(244)    = .true. ! outlier, lat = 75
346          longwatercaptag(246)    = .true. ! outlier, lat = 75
347          longwatercaptag(250)    = .true. ! outlier, lat = 75
348          longwatercaptag(252)    = .true. ! outlier, lat = 75
349          longwatercaptag(254)    = .true. ! outlier, lat = 75
350          longwatercaptag(256)    = .true. ! outlier, lat = 75
351          endif
352!--------------------------------------------------------------       
353
354
355           
356         else if (klon_glo .ne. 1) then
357       
358          print*,'No predefined ice location for this resolution :',
359     &           nbp_lon,nbp_lat-1
360          print*,'Please change icelocationmode in surfini.F'
361          print*,'Or add some new definitions ...'
362          call abort
363         
364         endif
365
366         do ig=1,klon_glo
367          if (longwatercaptag(ig)) watercaptag_glo(ig) = .true.
368         enddo
369
370
371        ELSE IF (icelocationmode .eq. 3) THEN
372     
373         print*,'Surfini: ice caps defined by lat and lon values'
374
375         do ig=1,klon_glo
376         
377c-------- Towards olympia planitia water caps -----------
378c--------------------------------------------------------
379
380          if ( ( ( lati_glo(ig)*180./pi .ge. 77.  ) .and. ! cap #2
381     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
382     .           ( long_glo(ig)*180./pi .ge. 110. ) .and.
383     .           ( long_glo(ig)*180./pi .le. 181. ) )
384     .         .or.
385
386     .         ( ( lati_glo(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
387     .           ( lati_glo(ig)*180./pi .le. 76.  ) .and.
388     .           ( long_glo(ig)*180./pi .ge. 150. ) .and.
389     .           ( long_glo(ig)*180./pi .le. 168. ) )
390     .         .or.
391     .         ( ( lati_glo(ig)*180./pi .ge. 77 ) .and. ! cap #5
392     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
393     .           ( long_glo(ig)*180./pi .ge. -150.) .and.
394     .           ( long_glo(ig)*180./pi .le. -110.) ) )
395     .         then
396             
397               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
398              !    watercaptag(ig)=.true.
399                  alternate = 1
400               else
401                  alternate = 0
402               endif !end if alternate = 0
403               
404          endif
405
406c----------- Opposite olympia planitia water cap --------
407c--------------------------------------------------------
408
409          if ( ( ( lati_glo(ig)*180./pi     .ge.  80 ) .and.
410     .         ( lati_glo(ig)*180./pi     .le.  84 ) )
411     .         .and.
412     .       ( ( long_glo(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
413     .         ( long_glo(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
414!     .     ( ( ( long_glo(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
415!     .         ( long_glo(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
416!     .       ( ( long_glo(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
417!     .         ( long_glo(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
418        !   watercaptag_glo(ig)=.true.
419          endif
420
421
422c -------------------- Central cap ----------------------
423c--------------------------------------------------------
424
425          if (abs(lati_glo(ig)*180./pi).gt.80)
426     .          watercaptag_glo(ig)=.true.
427           
428c--------------------------------------------------------
429c--------------------------------------------------------
430         end do ! of (klon_glo)
431
432
433        ELSE
434     
435         print*, 'In surfini.F, icelocationmode is ', icelocationmode
436         print*, 'It should be 1, 2 or 3.'
437         call abort
438
439        ENDIF ! of if (icelocation)
440       
441       
442        ! print caps locations - useful for plots too
443        print*,'surfini: latitude | longitude | ig'
444        do ig=1,klon_glo
445          dryness_glo(ig) = icedryness
446
447          if (watercaptag_glo(ig)) then
448             print*,'surfini: ice water cap', lati_glo(ig)*180./pi,
449     &              long_glo(ig)*180./pi, ig
450          endif
451        enddo
452       
453       endif !of if (is_master)
454       
455       ! Now scatter fields watercaptag and dryness from master to all
456       ! (is just a plain copy in serial mode)
457       call scatter(dryness_glo,dryness)
458       call scatter(watercaptag_glo,watercaptag)
459       
460! end of #else of #ifdef MESOSCALE
461       ENDIF ! (caps & water)
462#endif       
463
464c ===============================================================
465c      INITIAL ALBEDO
466c ===============================================================
467
468         write(*,*)"surfini: water frost thickness",
469     s     frost_albedo_threshold
470         write(*,*)"surfini: water ice albedo:", albedo_h2o_ice
471         write(*,*)"surfini: water ice TI:", inert_h2o_ice
472
473c        To start with : Initial albedo = observed dataset
474c        -------------------------------------------------
475         DO ig=1,ngrid
476              psolaralb(ig,1)=albedodat(ig)
477              psolaralb(ig,2)=albedodat(ig)
478         END DO
479         PRINT*,'surfini: minimum albedo without water caps',
480     &          minval(albedodat)
481         PRINT*,'surfini: maximum albedo without water caps',
482     &          maxval(albedodat)
483
484c        initial albedo if permanent H2O ice is present
485c        ------------------------------------------------
486         IF ((water) .and. (caps)) THEN
487           DO ig=1,ngrid
488            IF (watercaptag(ig)) THEN
489              psolaralb(ig,1) = albedo_h2o_ice
490              psolaralb(ig,2) = albedo_h2o_ice
491            ENDIF
492           END DO
493           PRINT*,'surfini: minimum albedo with water caps',
494     &            minval(albedodat)
495           PRINT*,'surfini: maximum albedo with water caps',
496     &            maxval(albedodat)
497
498         ENDIF
499
500c      changing initial albedo if CO2 ice is present
501c      -------------------------------------------
502
503       DO ig=1,ngrid
504         IF (piceco2(ig) .GT. 0.) THEN
505             IF(latitude(ig).LT. 0.) THEN
506                icap=2 ! Southern hemisphere
507             ELSE
508                icap=1 ! Northern hemisphere
509             ENDIF
510             psolaralb(ig,1) = albedice(icap)
511             psolaralb(ig,2) = albedice(icap)   
512         END IF
513       END DO
514
515c      changing initial albedo if water ice frost is present
516c      -------------------------------------------
517       IF (water) THEN
518          do iq=1,nqmx
519c          if there is frost and surface albedo is set to albedo_h2o_ice
520           if(noms(iq).eq."h2o_ice") then
521             do ig=1,ngrid
522             
523              if ((watercaptag(ig).eqv..false.)
524     &     .and. (qsurf(ig,iq).lt.-frost_albedo_threshold)) then
525              print*, ''
526              print*, '!!! PROBLEM in SURFINI !!!!'
527              print*, 'FOUND NEGATIVE SURFACE ICE VALUE WHERE
528     & WATERCAPTAG IS FALSE'
529              print*, ''
530              print*, 'ig,qsurf,threshold' ,
531     &         ig, qsurf(ig,iq), -frost_albedo_threshold
532              print*, ''
533              print*, '1) Check h2o_ice in startfi and ice
534     & distribution in surfini'
535              print*, '2) Use ini_h2osurf option in newstart'
536              print*, ''
537#ifndef MESOINI
538              CALL ABORT
539#else
540              qsurf(ig,iq) = 0.
541#endif
542            endif
543             
544              if ((piceco2(ig) .eq. 0.).and.
545     &          (qsurf(ig,iq).gt.frost_albedo_threshold)) then
546                     psolaralb(ig,1) = albedo_h2o_ice
547                     psolaralb(ig,2) = albedo_h2o_ice
548c                     PRINT*,'surfini.F frost',
549c     &                  latitude(ig)*180./pi, longitude(ig)*180./pi
550               endif
551              enddo
552           endif
553          end do
554          PRINT*,'surfini: minimum albedo with frost and co2',
555     &            minval(albedodat)
556          PRINT*,'surfini: maximum albedo with frost and co2',
557     &            maxval(albedodat)
558       ELSE
559         watercaptag(:) = .false.
560       END IF ! OF IF(water)
561         
562      RETURN
563      END
Note: See TracBrowser for help on using the repository browser.